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Analysis  of  a  System  to  Prevent  Helicopter 
Rotor  Blade-Airframe  Strikes 


B.  W.  McCormick1  and  R  G.  Melton2 
Abstract 

Rotor  blade-airframe  strikes  are  rare  but  they  do  occur.  Three  areas  of  the 
airframe  are  particularly  vulnerable:  the  tail  boom,  canopy  and,  in  the  case  of  the 
underslung,  teetoring  rotor,  the  rotor  shaft.  This  latter  case  is  known  as  mast  bumping. 
This  report  studies  a  system  to  prevent  a  helicopter  rotor  blade  from  striking  any  part 
of  the  airframe.  Essentially,  the  system  continuously  predicts  ahead  the  rotor  blade 
flapping  in  response  to  an  input  such  as  pilot  control  or  an  atmospheric  disturbance.  If 
a  blade  strike  is  predicted  to  occur  then  an  appropriate  feedback  control  is  applied  to 
alter  the  future  flapping.  The  prediction  is  then  begun  again  with  the  altered  control. 

In  the  actual  system,  an  enunciator  might  warn  the  pilot  at  the  time  that  he  is 
attempting  a  control  input  which  could  be  hazardous.  Two  somewhat  independent 
approaches  to  the  design  of  the  controller  are  taken.  One  of  the  programs  is  entirely 
numerical  in  its  approach.  The  other  utilizes  modem  control  theory  and  considers  the 
preliminary  aspects  of  implementing  the  controller  in  digital  hardware.  Both  methods 
indicate  the  feasibility  of  preventing  excessive  flapping,  although  the  question  of 
implementation  in  a  dedicated  microprocessor  is  not  fully  resolved. 

Nomenclature 

The  nomenclature  is  to  be  found  following  the  text.  It  includes  both  symbols  used 
in  the  text  as  well  as  those  used  in  the  computer  codes. 

Introduction 

The  purpose  of  this  study  is  to  determine  the  feasibility  of  a  control  system  which 
will  prevent  a  helicopter  rotor  blade  from  striking  any  part  of  the  airframe.  Essentially, 
the  idea  of  the  system  is  to  continually  predict  ahead  how  the  rotor  blade  will  flap  in 
response  to  an  input  to  the  rotor  such  as  pilot  control  or  an  atmospheric  disturbance.  If 
a  blade  strike  is  predicted  to  occur  then  an  appropriate  feedback  control  is  applied  to 
alter  the  future  flapping.  The  prediction  is  then  begun  again  with  the  altered  control. 

In  the  actual  system,  an  enunciator  might  warn  the  pilot  at  the  time  that  he  is 
attempting  a  control  input  which  could  be  hazardous. 
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Rotor  blade-airframe  strikes  have  occurred  on  several  helicopter  configurations. 

Three  areas  of  the  airframe  are  particular  vulnerable,  the  tail  boom,  canopy  and,  in  the 
case  of  the  underslung,  teetoring  rotor,  the  rotor  shaft.  This  latter  case,  known  as 
"mast  bumping",  was  the  primary  motivation  for  this  study  but  the  results  should  be 
generally  applicable  to  other  helicopter  configurations. 

Two  somewhat  independent  approaches  to  the  design  of  the  controller  have  been 
taken.  One  of  the  programs  is  entirely  numerical  in  its  approach.  The  other  utilizes 
modern  control  theory.  Both  approaches  require  that  the  following  questions  be 
answered. 

1.  How  many  revolutions  of  the  rotor  must  the  flapping  be  predicted  ahead? 

2.  What  type  of  feedback  control  is  required? 

3.  How  quickly,  in  terms  of  a  rotor  revolution,  must  the  microprocessor  predict 

the  future  flapping? 

The  answers  to  these  questions  must  depend,  in  part,  upon  the  operating  state  of 
the  rotor,  the  point  around  the  azimuth  where  the  strike  is  predicted  to  occur,  and  the 
severity  of  the  predicted  strike. 

The  material  to  be  developed  here  attempts  to  answer  the  above  questions  but  with 
certain  limitations.  It  proved  to  be  a  more  difficult  task  than  had  been  anticipated  so 
that,  with  the  allotted  funding,  it  has  not  been  possible  to  include,  as  yet,  the  dynamics 
of  the  airframe,  nor  to  determine,  with  any  certainty,  whether  or  not  a  microprocessor 
can  provide  the  speed  which  is  necessary.  Based  on  informal  discussions  with  persons 
knowledgeable  in  the  design  of  microprocessors,  it  is  felt  at  this  time  that  the  necessary 
speed  can  be  achieved.  Also,  it  may  be  possible,  within  the  accuracy  required,  to  replace 
some  of  the  numerical  integrations  with  approximate  closed-form  expressions,  thus 
increasing  the  speed  of  the  calculations. 

This  report  begins  by  presenting  the  development  of  a  non-linear,  numerical  program 
to  predict  the  flapping  of  a  rotor  including  retreating  blade  stall  and  reversed  flow. 

This  program  is  utilized  in  both  the  numerical  study  and  in  the  one  utilizing  modern 
control  theory.  Next,  the  logic  for  the  numerical  controller  is  presented  together  with 
some  results  which  were  obtained  using  the  AH-1J  helicopter  as  an  example.  This  is 
followed  by  the  analytical  developments  based  on  modern  control  theory  and  some  results 
of  that  analysis,  again  using  the  AH-1J  as  an  example.  Finally,  some  conclusions  and 
recommendations  are  made,  the  principal  one  being  that  the  scheme  for  preventing  blade 
strikes  appears  to  be  feasible  and  should  be  pursued  further. 

Description  of  Program  to  Predict  Rotor  Blade  Flapping 

A  program  to  predict  blade  flapping  was  written  specifically  for  this  effort  for  two 
reasons.  First,  it  was  uncertain  that  a  classical  approach  which,  at  any  azimuth  position, 
obtains  the  integrated  blade  lift  and  hub  moment  in  closed  form,  would  be  adequate.  The 
classical  approach  is  limited  to  first  harmonic  flapping  and  contains  certain  small  angle 
assumptions  which  may  be  significant.  The  classical  approach  also  neglects  reverse  flow 
and  retreating  blade  stall.  Secondly,  at  the  other  extreme,  it  was  felt  that  the 
computational  time  required  by  existing  elaborate  codes  which  predict  rotor  flapping,  such 
as  C-81,  would  be  prohibitive  for  the  proposed  control  system.  Thus  a  compromise 
between  these  two  extremes  was  taken. 


3 


The  subroutine  which  was  written  begins  with  the  rotor  state  at  a  particular  instant 
and  azimuth  angle,  1>,  and  integrates  the  thrust  and  torque  over  the  radius.  A  uniform 
downwash  is  assumed  together  with  a  rotor  blade  which  is  rigid  but  semi -articulated  with 
only  a  flapping  degree  of  freedom.  A  better  model  of  the  downwash,  like  a  triangular 
variation  or  a  prescribed  wake,  could  be  incorporated  into  this  model  but  consideration 
of  a  flexible  blade  would  require  extensive  modification.  At  every  azimuth  position,  the 
blade  loading  is  numerically  integrated  along  the  span  to  obtain  the  instantaneous  lift 
and  hub  moment.  No  small  angle  assumptions  are  made  with  regard  to  the  angle  of 

attack  of  the  blade  sections.  Reverse  flow  and  stall  are  accounted  for  by  the  use  of  a 
table  lookup  to  obtain  the  airfoil  lift  and  drag  coefficients  for  angles  of  attack  from 
zero  to  360  degrees. 

For  completeness,  the  analytical  basis  for  the  flapping  program  will  be  presented  in 
detail.  Consider  figure  1  which  is  a  left  side  view  of  a  rotor  with  the  disc  plane  at  an 
angle  of  attack.  The  disc  plane  is  sometimes  referred  to  as  the  shaft  plane  and  is  the 
plane  normal  to  the  shaft.  From  this  figure,  it  is  seen  that  the  freestream  velocity,  V, 
can  be  split  into  two  components,  one  normal  to  the  disc  plane  and  the  other  lying  in 
the  plane. 


Now  consider  figure  2.  This  figure  is  a  top  view  of  the  disc  with  the  blade  at  an 
azimuth  angle,  rp,  measured  clockwise  from  the  downstream  position.  It  is  seen  that  the 
in-plane  velocity  can  be  further  divided  into  two  components,  one  parallel  and  the  other 
normal  to  the  blade.  Also  shown  in  this  figure  is  the  linear  velocity  component  directed 
normal  to  the  blade  at  a  radius  of  r  which  results  from  the  angular  velocity  of  the  rotor. 


Figure  3  is  a  view  in  the  plane  defined  by  the  blade  and  the  shaft  axis.  The  blade 
is  shown  with  a  flapping  angle,  P,  a  flapping  velocity,  dP/dt  and  an  angular 
acceleration  about  the  flapping  axis.  At  a  radius  of  r,  if  the  blade  is  flapping  up,  as 
shown  relative  to  the  blade,  a  downward  velocity  results  from  the  upward  flapping. 
Further,  the  previous  component  of  the  velocity  parallel  to  the  rotor  has,  itself,  a 
component  directed  up  normal  to  the  blade. 

Figure  4  shows  the  blade  section  at  the  radius,  r,  at  a  pitch  angle  6  relative  to 
the  disc  plane.  The  net  velocity  up,  normal  to  the  disc  plane,  is  given  by. 


Where: 


V#  «  V  sin  -  w  -  (r-c)/?  -  0V  cos  «*;  cos  i> 

w  «  rotor  downwash  velocity 
r  ■  radial  distance  of  blade  section  from  shaft  axis 
e  m  radial  distance  of  flapping  hinge  from  shaft  axis 
^  -  azimuth  angle 

V  »  velocity  at  which  rotor  is  advancing 
p  ■  flapping  angle  between  rotor  blade  and  disc  plane 
a;  *  disc  plane  angle  of  attack  relative  to  V 


The  net  velocity  in  the  disc  plane  can  be  obtained  from. 


Where: 


Vt  -  wr  +  Vcos  sin 
w  «  angular  velocity  of  rotor 


0) 


(2) 


Thus,  from  figure  4  and  equations  (1)  and  (2),  the  angle  of  attack  of  the  blade 
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section  is  given  by. 


Where: 

4  -  tan  -\Vt/Vu) 


a  =  0  +  ^ 


(3) 


One  must  be  careful  in  programming  the  above  to  remember  that  the  angle  4  can 
lie  in  any  quadrant  and  can  be  of  a  magnitude  such  that  the  angle  of  attack,  of 
the  section  can  vary  from  0  to  360  degrees.  In  order  to  obtain  the  correct  section 
Ci  and  Cd,  and  to  correctly  resolve  the  lift  and  drag  forces  into  the  thrust  and 
torque  directions,  one  should  check  the  sign  of  both  Vt  and  Vu.  Figure  5  illustrates 
the  possible  flow  directions  and  the  corresponding  lift  and  drag  vectors.  In  the 
subroutine  SUBFLAP,  which  is  appended  to  this  report  and  which  will  be  discussed  in 
more  detail  later,  in  lieu  of  using  equation  (3)  to  calculate  the  angle  of  attack,  the  angle 
4  is  first  calculated  simply  as: 


*  -  tan  w 

The  angle  of  attack  is  then  found  by  the  four  possible  combinations  shown  in  figure  5. 
Specifically, 


V, 

> 

0 

and 

vu 

> 

0 

oc  *  ©  +  ^ 

(5a) 

V, 

> 

0 

and 

Vu 

< 

0 

a  b  0  -  <f> 

(5b) 

vt 

< 

0 

and 

Vu 

> 

0 

OC-5T  +  0- 

4 

(5c) 

< 

0 

and 

Vu 

> 

0 

a  *  x  -  0  - 

4 

(5d) 

In  the  normal  case,  the  lift  coefficient,  Cj,  adds  to  the  rotor  thrust  and  to 
the  torque  while  the  drag  coefficient,  Cd,  adds  to  the  torque  and  subtracts  from  the 
thrust.  However,  in  the  reverse  flow  region,  or  anywhere  along  the  blade  where  the 
angle  of  attack  is  greater  than  90  degrees,  this  is  no  longer  the  case.  Thus  factors 
multiplying  C j  and  Cd  are  set  equal  to  1  or  -1  depending  upon  whether  or  not 

the  resultant  flow  is  impinging  on  the  leading  edge  or  trailing  edge  and  from  above  or 
below. 

Data  for  airfoils  over  an  alpha  range  from  0  to  360  degrees  is  difficult  to  find  so 
that,  for  this  study,  the  data  for  the  NACA  0012  airfoil  from  0  to  180  degrees  was  used. 
This  is  the  same  data  used  in  C-81.  The  logic  for  determining  C^,  Cd  and  the 
factors  for  resolving  the  lift  is  found,  and  identified,  in  the  appended  subroutine. 

Knowing  Cj  and  Cd  from  a  table  lookup  and  the  angles  above,  the 
derivatives  of  the  rotor  lift  and  drag  can  be  found  from, 

^■jpVe  cC i  (6a) 

®  ■  j  p  Ve  c  (6b) 

where  V(  is  the  resultant  velocity  shown  in  figure  4  and  given  by, 

V,  -  (  Vu  +  Vt  )  1/2 

The  radial  derivative  of  the  blade  thrust  is  then  found  from. 


(7) 
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dT  dL  I  vt  I  dD  I  vu  I 
<JF  =  dr  Vs  +  dr  V, 


(8) 


The  derivative  of  the  moment  about  the  flapping  hinge  is  found  by  multiplying  the  above 
by  the  distance,  r-e,  from  the  flapping  hinge  to  the  blade  element.  Thus, 


w-  <r  - S  <»> 

Using  a  simple  trapezoidal  rule,  equations  (8)  and  (9)  are  integrated  along  the  blade 
from  e  to  R  to  obtain  the  total  instantaneous  thrust  and  hinge  moment.  Knowing  the 
blade  moment  of  inertia  and  the  moment  of  the  blade  weight  about  the  flapping  hinge, 
the  angular  acceleration  about  the  flapping  hinge  is  then  calculated  from. 


'rjj?  -  M  -  IF  P  w*  -  Mw 


(10) 


The  angular  velocity  and  the  flapping  angle,  p,  at  the  end  of  the  time  increment,  A  t, 
are  then  obtained  from. 


dp 


dp 


d2p 


g  ft  ♦  A  t)  -  £  (t)  ♦  At 

f  (•  *  *>  ■  *<>  ♦  Sr*  ♦  0 


(ID 

(12) 


This  is  essentially  the  end  of  the  flapping  subroutine.  The  new  state  of  the  rotor 
defined  by  (11)  and  (12)  and  the  instantaneous  thrust  from  (8)  is  returned  to  the  main 
program  to  be  integrated  with  respect  to  the  azimuth  angle,  V>,  and  to  be  used  in  the 
logic  of  the  controller. 


Numerical  Controller 


The  action  of  the  controller  is  shown  schematically  in  figure  6.  Here,  the  flapping 
angle,  P,  is  shown  as  a  function  of  time.  Suppose,  for  example,  that  the  controller  is 
activated  at  the  time  corresponding  to  point  A.  At  that  time  it  begins  to  predict  the 
flapping  of  the  rotor  ahead  for  a  specified  number  of  revolutions  based  on  the  state  of 
the  rotor  at  the  time  and  on  a  linear  extrapolation  of  the  control  input.  If  the  rotor 
flapping  is  predicted  to  be  within  limits,  the  controller  returns  to  the  rotor  after  a  time, 
r,  which  is  the  time  required  to  perform  the  calculations.  During  this  time  the  rotor 
blade  has  moved  from  A  to  point  B.  The  process  is  then  repeated.  As  illustrated,  the 
flapping  at  point  A  was  predicted  to  be  within  prescribed  limits  for  three  revolutions 
ahead.  Thus  the  controller  simply  returns  after  the  time  r  to  sense  a  new  rotor  state 
in  order  to  perform  another  prediction. 

Suppose,  at  some  later  time,  point  C,  the  controller  begins  a  prediction  during 
which,  because  of  a  control  input,  the  rotor  flapping  is  predicted  to  exceed  a  limiting 
value.  Upon  reaching  this  limit  the  prediction  immediately  stops  and  the  controller 
returns  to  the  rotor  at  some  fraction  of  r,  point  D.  At  this  instant  it  commands  an 
incremental  step  input  of  corrective  control  to  prevent  the  excessive  flapping  which  was 
predicted.  It  then  begins  a  new  prediction  with  the  incremental  control.  If  excessive 
flapping  is  not  predicted  then  an  incremental  step  of  corrective  control  is  removed. 
Thus,  the  control  actuators  commanded  by  the  controller  are  generally  inactive  except 
for  those  rare  occasions  when  excessive  control  is  applied  by  the  pilot  or  when  external 
disturbances  are  sufficient  to  cause  excessive  flapping. 
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The  logic  for  the  above  is  shown  in  figure  7.  This  logic  serves  several  purposes. 

First,  it  calls  for  the  subroutine  FLAP  in  order  to  model  the  operating  rotor.  Secondly 
it  calls  for  this  same  subroutine  to  model  the  action  of  predicting  the  future  flapping  of 
the  rotor  and  finally,  it  checks  the  future  flapping  against  prescribed  limits  and  provides 
the  necessary  feedback  control  to  alter  the  predicted  future  flapping.  The  subroutine 
DNWSH  provides  the  downwash  velocity  as  a  function  of  blade  azimuth  position  and  rotor 
thrust.  To  date,  as  stated  previously,  only  a  simple  uniform  downwash  model  has  been 
used.  The  subroutine  CNTRL  provides  the  corrective  control  which  is  a  function  of  how 
excessive  the  flapping  is  predicted  to  be  and  the  azimuth  position  at  which  it  is 
predicted  to  occur. 

Figure  7  will  now  be  explained  in  some  detail.  The  algorithm  begins  by  inputting 
the  rotor  state,  operating  conditions  and  parameters  defining  the  controller.  The  state 
of  the  rotor;  ie,  the  advance  velocity,  trim  angle  of  attack  and  control  angles  are  read 
from  data  files  with  the  angles  being  determined  from  the  static  trim  program  described 
in  reference  (1).  Parameters  governing  the  rotor  flapping,  namely  the  rotor  geometry 
and  inertia  properties,  are  also  read  in  from  data  files.  The  variables  specific  to  a 
particular  numerical  calculation  are  read  in  from  the  keyboard  by  the  operator. 
Specifically,  in  order,  these  are: 

1.  Identifying  case  number 

2.  Maximum  flapping  angle  to  be  allowed,  BETLIM 

3.  Increments  to  cyclic  pitches  if  allowable  flapping  is  predicted  to  be 
exceeded,  DELFB1  and  DELFB2 

4.  Maximum  feedback  on  cyclic  controls,  FB1LIM  and  FB2LIM 

5.  Number  of  rotations  before  control  input  (to  disturb  static  trim) 
is  applied,  N 

6.  Number  of  rotations  that  flapping  is  to  be  predicted  ahead  at  a  given 
instant,  NPRED 

7.  Fraction  of  a  revolution  required  to  accomplish  the  above  prediction,  FPRED 

8.  Rate  of  linear  increase  of  cyclic  control  input  to  disturb  system,  CRATE1 
and  CRATE2 

9.  Maximum  incremental  values  for  the  cyclic  control  inputs,  THE1LM 
and  THE2LM 

1  0  .  Length  of  time  for  cyclic  control  inputs  to  be  applied,  TOFF 

1  1  .  Number  of  numerical  integrations  between  printouts,  PRT 

1  2  .  Maximum  number  of  revolutions  for  run,  NMAX 

Following  the  input,  the  time,  azimuth  angle  and  quantities  to  be  integrated  are 
initialized  to  either  zero  or  to  their  initial  values.  The  switch  KNTRL  is  set  to  1  and 
the  calculations  begin.  The  circled  nodes  in  Figure  7  numbered  1  through  21  refer  to 
corresponding  statement  numbers  in  the  program  FLAP.  The  value  of  KNTRL  remains  at 
1  until  the  time  is  reached  for  the  control  input  to  begin.  At  any  instant  of  time  and 
azimuth  angle,  i>,  the  subroutine  SUBFLAP  integrates  the  blade  thrust  with  respect  to 
radius  in  order  to  obtain  the  derivative  of  the  total  rotor  thrust  with  respect  to  i>. 

Following  statement  5,  V>  is  checked  to  see  if  it  exceeds  the  value  of  2*  or  the 
value  at  which  control  input  begins.  If  so,  KNTRL  is  set  to  2  and,  after  saving  the 
current  state  variables  as  initial  values  for  later  use,  the  prediction  of  the  future  rotor 
flapping  begins. 
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The  prediction  continues  over  NPRED  revolutions  or  until  the  predicted  flapping 
exceeds  BETLIM.  During  all  of  the  calculations,  after  every  revolution,  the  thrust  is 
averaged  and  the  subroutine  DWNWSH  called  to  update  the  average  downwash  velocity. 
The  time,  TC,  required  for  the  predictions  is  calculated  following  statement  14, 
correspondent  to  the  time  at  which  BETLIM  was  exceeded,  or  simply  equated  to  the  time 
required  for  NPRED  revolutions  if  BETLIM  was  not  exceeded.  The  program  then  sets 
K.NTRL  equal  to  3  and  returns  to  the  initial  conditions  saved  when  the  prediction  started. 
Calculations  of  the  flapping  are  then  repeated  until  the  time,  TC,  is  reached.  At  this 
point,  the  parameter  FB  is  increased  by  1  which  adds  an  increment  of  feedback  control. 

If  it  is  predicted  that  BETLIM  will  not  be  exceeded,  then  FB  is  reduced  by  1.  However, 
the  subroutine  SUBCNTRL  does  not  allow  FB  to  be  less  than  zero.  At  this  time,  and 
with  the  feedback,  KNTRL  is  set  again  to  zero  and  the  prediction  begins  again.  This 
alternating  process  of  predicting  the  future  flapping  and  calculating  the  actual  flapping 
with  feedback  continues  until  NMAX  revolution  are  reached. 

Results  from  the  Numerical  Program 

The  AH-1J  helicopter  was  used  as  an  example  to  demonstrate  the  programs 
developed  here  with  numerical  values  of  the  parameters  for  this  helicopter  and  its  rotor 
system  being  obtained  from  reference  (2). 

Comparison  with  Predictions  from  Classical  Theory 

To  begin,  a  comparison  was  made  between  predictions  of  steady  flapping  based  on 
the  numerical  code  with  those  based  on  classical  theory,  (see  reference  3)  The  typical 
results  shown  in  Figures  8  and  9  appear  to  confirm  the  code  and  offer  the  promise  of 
being  able  to  use  the  linearized  approximations  contained  in  the  classical  formulation  to 
speed  up  the  numerical  controller.  For  these  particular  operating  states,  the  amplitude 
of  the  flapping,  as  predicted  by  the  classical  theory,  agrees  almost  exactly  with  the 
predictions  from  the  numerical  model.  These  particular  graphs  are  for  an  AH-1J 
operating  at  61  knots  at  3075  feet  at  a  gross  weight  of  9500  lbs.  In  figure  8,  the 
controls  are  held  constant.  In  Figure  9,  the  lateral  cyclic  is  increased  rapidly  by  10 
degrees  after  two  revolutions  and  then  held  constant. 

Rotor  Response 

The  response  of  the  rotor  to  a  lateral  control  input  is  presented  in  Figure  10. 

Here,  the  helicopter  is  trimmed  at  80  kts  at  standard  sea  level  conditions  using  the  static 
trim  program  in  reference  1.  The  rotor  revolutions  are  measured  starting  with  the  blade 
at  an  azimuth  angle  of  zero.  Since  this  latter  program,  based  on  classical,  linearized 
theory  differs  slightly  from  the  numerical  predictions  of  flapping,  it  takes  approximately 
two  revolutions  for  the  numerical  results  to  become  steady.  At  3.0  and  at  3.5 
revolutions,  a  step  increment  in  the  lateral  cyclic  control  input  of  5  degrees  is  applied. 

It  can  be  seen  that  the  flapping  responds  rapidly  to  the  control  input  and  reaches  a 
steady  state  in  approximately  one  and  one-half  revolutions  or  less.  When  the  lateral 
cyclic  control  is  applied  at  an  azimuth  angle  of  zero  degrees  (3  revolutions),  the  lateral 
flapping  is  seen  to  increase  by  approximately  two  degrees  only  a  quarter  of  a  revolution 
later  at  a  psi  value  of  90  degrees,  and  by  three  and  one-half  degrees  another  half  of  a 
revolution  later  at  a  psi  of  270  degrees. 

The  effect  of  delaying  a  correction  to  a  disturbance  is  illustrated  by  Figures  11  and 
12.  In  both  figures,  a  lateral  cyclic  control  of  5  degrees  is  applied  at  the  end  of  3 
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revolutions  (1080  degrees).  Then,  in  Figure  11  a  correction  of  -5  degrees  is  applied  0.4 
of  a  revolution  later  whereas,  in  Figure  12,  the  correction  is  not  applied  until  one 
revolution  later.  This  same  case  was  run  for  correction  delays  of  0.1,  0.2,  and  0.3  of  a 
revolution  with  results  which  were  nearly  identical  to  figure  11.  It  can  be  seen  that  the 
quicker  response  in  figure  11  does  not  appreciably  reduce  the  peak  flapping  immediately 
after  the  input  but  the  next  peak  is  reduced  by  approximately  2  degrees  in  comparison  to 
the  slower  response  of  Figure  12. 


Effect  of  Feedback 

Figures  8  through  12  typify  the  kinds  of  evaluations  which  were  done  prior  to 
constructing  the  numerical  controller  previously  discussed.  After  gaining  confidence  in 
the  numerical  flapping  program  and  getting  a  feel  for  the  feedback  capability  which 
would  be  needed,  the  final  program  was  developed  and  parametric  studies  performed.  To 
date,  these  have  only  been  done  for  the  AH-1J  at  61  knots,  3075  ft.  and  9500  lbs.  ft. 
Some  typical  results  illustrating  the  effects  of  the  various  parameters  are  shown  in 
Figures  13  through  17.  Table  1  lists  the  parameters  used  for  each  case  noted  on  the 
figures.  The  parameters  which  are  varied  in  these  figures  include: 

a.  The  rate  at  which  the  controls  can  be  moved. 

b.  The  feedback  increment  to  the  controls  each  time  excessive  flapping  is 
predicted. 

c.  The  time  required  to  predict  the  future  flapping  as  a  fraction  of  a  rotor 
revolution. 

d.  The  limit  on  the  authority  of  the  feedback  controls. 

e.  The  number  of  revolutions  ahead  through  which  the  future  flapping  is 

predicted. 

Figure  13  illustrates  the  effect  of  the  rate  of  control  movement.  To  put  all  rates 
and  times  in  perspective,  it  is  noted  for  the  AH-1J  that  it  takes  0.19  seconds  per 

revolution  or  5.2  revolutions  per  second.  Thus,  for  example,  at  the  rate  shown  in  Figure 

13  of  100  degs/sec,  the  controls  will  move  19.2  degrees  in  one  revolution.  From  this 

figure,  it  would  appear  that  the  application  rate  of  the  controls  is  not  a  predominant 

effect.  It  is  emphasized  that  this  is  the  rate  at  which  the  controls  are  being  moved  to 
disturb  the  flapping  and  is  not  related  to  the  rate  of  control  feedback.  The  latter  effect 
is  accomplished  in  steps  and  is  discussed  later.  It  should  be  noted  that  the  flapping 
limit  for  these  calculations  and  for  all  of  the  results  to  follow  was  set  at  8  degrees. 

Also,  the  disturbance  for  all  of  the  cases  was  produced  by  applying  10  degrees  of  lateral 

cyclic  at  the  end  of  two  revolutions.  In  the  case  of  Figure  13  it  appears  that  the 

combination  of  looking  ahead  only  two  revolutions,  feedback  increments  of  2  degrees,  and 
an  FPRED  value  of  0.2  results  in  a  flapping  which  will  exceed  the  flapping  limit  slightly. 
The  err  t  of  these  parameters  are  discussed  in  the  following  paragraphs. 

Ti  feedback  controls  are  applied  as  an  accumulation  of  step  inputs.  A  step  is 
appiiea  each  time  the  controller  predicts  excessive  flapping  and  is  removed  each  time  the 
controller  "edicts  that  the  flapping  will  remain  within  limits.  The  effect  of  the 
fee  'back  •  tep  size  is  illustrated  in  Figure  14.  Here  again,  the  influence  of  the  step  size 
is  not  a  dominant  parameter,  at  least  for  this  particular  operating  state.  However,  it 
does  appear  as  if  a  feedback  increment  of  at  least  4  degrees  is  needed  to  stay  within  the 
flapping  limit. 
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The  parameter,  FPRED,  is  the  fraction  of  a  revolution  required  to  predict  the  future 
flapping  of  the  rotor  over  a  given  number  of  revolutions.  Figure  15  illustrates  the  effect 
of  varying  this  parameter  over  a  range  from  0.1  to  0.4  where  the  flapping  is  being 
predicted  three  revolutions  ahead.  The  feedback  step  size  for  this  figure  is  4  degrees. 

The  results  shown  here  are  somewhat  puzzling  and  tend  to  contradict  intuition.  One 
would  think  that  the  lower  FPRED  the  better.  However,  the  graph  shows  a  higher 
flapping  for  lower  values  of  FPRED  at  around  six  revolutions.  (although  initially 
following  the  disturbance  at  around  three  revolutions  the  results  are  opposite)  It  may  be 
that,  what  amounts  to  a  high  gain  in  the  feedback,  may  be  causing  some  type  of  an 
instability  in  the  flapping. 

Figure  16  shows  the  effect  of  limiting  the  total  feedback  control  angle.  From  this 

figure  it  appears  as  if  the  authority  of  the  feedback  cannot  be  limited  very  much.  For 

this  particular  case,  as  a  result  of  the  10  degrees  of  lateral  cyclic  input,  a  total  feedback 
of  approximately  7  degrees  must  be  allowed  in  order  to  stay  within  the  prescribed 
flapping  limit. 

The  effect  of  looking  ahead  1,  2,  and  3  revolutions  in  predicting  the  future  flapping 
is  shown  in  Figure  17.  As  one  might  expect,  the  differences  in  the  results  show  up  only 
in  the  first  couple  of  peaks  following  the  control  input.  Looking  ahead  only  one 
revolution  results  in  a  flapping  peak  at  about  2.8  revolutions  which  is  approximately  one 
degree  higher  than  for  the  other  two  cases.  It  appears,  from  these  and  other  cases,  as 

if  a  value  of  NPRED  of  2  or  higher  might  be  satisfactory. 
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Nomenclature  for  the  Numerical  Controller  Studies 

The  symbols  given  here  are  specific  to  the  foregoing  material  since  the  study  of  the 
application  of  modern  control  theory  was  essentially  independent  of  that  for  the 
numerical  controller.  A  separate  table  of  nomenclature  giving  the  additional  symbols  for 
the  analysis  using  modern  control  theory  can  be  found  at  the  end  of  that  material. 

Nomenclature  for  Text 

SYMBOL  DEFINITION 


« 

P 

At 

e 

<t> 

e 

V- 

w 

c 

Cd 

C/ 

D 

L 

M 

Mw 

r 

t 

T 

V 


Blade  section  angle  of  attack 
Rotor  disc  plane  angle  of  attack 
Blade  flapping  angle  measured  from  disc  plane 
Increment  in  time 

Distance  from  shaft  axis  to  flapping  hinge 

Angle  of  resultant  velocity 

Blade  section  pitch  angle 

Blade  azimuth  angle 

Angular  velocity  of  rotor 

Blade  section  chord 

Blade  section  drag  coefficient 

Blade  section  lift  coefficient 

Blade  section  drag 

Blade  moment  of  inertia  about  flapping  axis 
Blade  section  lift 

Aerodynamic  moment  about  flapping  axis 
Blade  weight  moment  about  flapping  axis 
radial  distance  from  shaft  axis 
time 

Blade  thrust 
Advance  velocity 


Vc  Resultant  velocity  at  blade  section 

Vt  Tangential  velocity  at  blade  section 

Vu  Upward  velocity  at  blade  section 


Variable  Name 

A1 

AI 

ALPHA 

ALPHAB 

ALPHAD  (I) 

ALPHAL  (I) 

AREA 

B 

B1 

BETA 

BETAO 

BETAD 

BETADG 

BETADI 

BETAI 

BETLIM 

BI 

C 

CO 

CALPHA 

CASE 

CBAR 

CD 

CDFAC 
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Nomenclature  for  Programs 
Definition 

longitudinal  flapping 

temporary  dummy  variable  in  interpolation  of  airfoil  table 

disc  plane  angle  of  attack 

blade  section  angle  of  attack 

section  drag  coefficient  tabulated  vs.  this  angle 

section  lift  coefficient  tabulated  vs.  this  angle 

disc  area  of  rotor 

number  of  blades 

lateral  flapping 

blade  flapping  angle  relative  to  disc  (shaft  axis)  plane 
coning  angle 

first  derivative  with  respect  to  time 

beta  in  degrees  for  printout 

initial  value  of  BETAD  when  starting  calculation 

initial  value  of  BETA  when  starting  prediction  of  flapping 

limiting  value  which  BETA  is  not  to  exceed 

temporary  dummy  variable  in  interpolation  of  airfoil  table 

blade  section  chord 

blade  section  root  chord  (X  +  0) 

cosine  of  alpha 

identifying  number 

mean  chord  of  linearly  tapered  blade 

section  drag  coefficient 

factor  resolving  drag  in  proper  direction  for  extreme  alpha 
tabulated  section  drag  coefficient  as  function  of  ALPHAD  (I) 


CDI  (I) 


Cl 

CL 

CLFAC 

CLI 

CPSI 

CRATE  1 

CRATE2 

CT 

D 

DDDR 

DELBT1 

DELBT2 

DELFB1 

DELFB2 

DELPSI 

DELR 

DELT 

DELT1 

DELT2 

DELTH1 

DELTH2 

DELTPR 

DELX 

DI 

DLDR 

DMDR 


temporary  dummy  variable  in  interpolation  of  airfoil  table 
section  lift  coefficient 

factor  resolving  lift  in  proper  direction  for  extreme  alpha 
tabulated  section  lift  coefficient  as  function  of  ALPHAD  (I) 
cosine  of  PSI 

rate  of  increase  of  lateral  cyclic,  degs/sec 
rate  of  increase  of  longitudinal  cyclic,  degs/sec 
rotor  blade  tip  chord 
rotor  diameter 

derivative  of  section  drag  with  respect  to  radius 
amount  by  which  BETA  exceeds  BETLIM 
temporary  vale  of  delbtl  in  order  to  determine  max  value 
increment  to  THETA  1  when  BETA  is  predicted  to  exceed  limit 
increment  to  THETA2  when  beta  is  predicted  to  exceed  limit 
increment  in  azimuth  angle  for  numerical  integration 
increment  in  radius  for  numerical  integration 
increment  in  time  for  numerical  integration 
differential  thrust  as  a  function  of  radius  for  a  given  PSI 
same  as,  and  averaged  with,  DELT1  to  integrate  for  thrust 
derivative  of  total  blade  thrust  with  azimuth  angle,  PSI 
same  as,  and  averaged  with,  DELTH1  to  integrate  for  TAVG 
increment  in  time  for  printout 

increment  in  dimensionless  radius  for  numerical  integration 
temporary  dummy  variable  in  interpolation  of  airfoil  table 
derivative  of  blade  lift  with  radius 

radial  derivative  of  blade  aerodynamic  moment  8t  flap  hinge 


DTDR 


radial  derivative  of  blade  thrust 
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I 


» 


DTR 

DWNWSH 

EPS 

EPSR 

FB 

FPRED 

KBETA 

KNTRL 

LAMDA 

Ml 

M2 

MAERO 

MIF 

MU 

MW 

N 

NMAX 

NPRED 

OMEGA 

PERIOD 

PI 

PRT 

PSI 

PSIDG 

PSII 

PSIII 


factor  to  convert  from  degrees  to  radians 

subroutine  to  calculate  average  downwash 

dimensionless  distance  of  flapping  hinge  from  rotor  axis 

distance  of  flapping  hinge  from  rotor  axis 

switch  to  add  or  subtract  feedback  correction  to  theta 

fraction  of  revolution  required  to  predict  NPRED  ahead 

rotor  pitch-flap  coupling  (delta-3) 

logic  switch,  see  Figure  7 

inflow  ratio  equals  net  flow  up  divided  by  tip  speed 

radial  derivative  of  aerodynamic  moment  about  flapping  hinge 

same  as,  and  averaged  with.  Ml  to  integrate  for  moment 

moment  about  flapping  axis  produced  by  aerodynamic  forces 

blade  mass  moment  of  inertia  about  flapping  axis 

ration  of  forward  speed  to  tip  speed,  V/VT 

blade  weight  moment  about  flapping  hinge 

number  of  rotations  at  which  control  initiated 

maximum  number  of  revolutions  for  run 

number  of  rotations  to  predict  flapping  ahead 

rotational  velocity  of  rotor,  radians/sec 

time  for  one  rotor  rotation 

the  usual  constant,  3.14159 

number  of  calculations  between  printouts 

azimuth  angle 

azimuth  angle,  degrees 

initial  value  of  PSI  (see  Figure  7) 

initial  value  of  PSI  (see  Figure  7) 

azimuth  angle  at  which  control  initiated 


PSILIM 


PSIMAX 

R 

RHO 

SALPHA 

SIGMA 

SPSI 

SUBCNTRL 

SUBFLAP 

TAU 

TAVG 

TC 

TCALC 

TCON 

TCONI 

TH1DG 

TH2DG 

THE1I 

THE1LM 

THE2I 

THE2LM 

THETA 

THETAO 

THETA 1 

THETA2 

THETAT 

THRUST 

TIME 


azimuth  angle  corresponding  to  NMAX 
rotor  radius 
air  mass  density 
sin  of  ALPHA 

rotor  section  solidity  equals  B*C/PI/R 
sin  of  PSI 

subroutine  provides  a  control  input  as  a  function  of  TCON 

subroutine  integrates  over  R  at  PSI  for  flapping  acceleration 

parameter  equals  MW/MIF/OMEGA**2 

average  rotor  thrust  in  one  revolution 

total  elapsed  time  for  predictin  of  future  flapping 

integral  of  the  rotor  thrust  with  PSI  for  one  revolution 

elapsed  time  from  when  the  control  is  first  applied 

time  when  control  is  initiated 

THETA  1  in  degrees  for  printout 

THETA2  in  degrees  for  printout 

initial  value  of  THETA  I 

maximum  incremental  value  for  lateral  control 

initial  value  of  THETA2 

maximum  incremental  value  for  longitudinal  control 

blade  section  pitch  angle 

initial  trim  collective  pitch 

initial  trim  lateral  cyclic  pitch 

INITIAL  trim  longitudinal  cyclic  pitch 

total  blade  twist 

instantaneous  total  blade  thrust  at  a  given  value  of  PSI 
elapsed  time  from  start  of  run 
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! 

I 


TIMEI 

TOFF 

TPRINT 

TWOPI 

V 

VR 

VT 

VTHETA 

VU 

W 

WI 

X 

XH 


time  at  which  prediction  started  of  future  flapping 
length  of  time  for  control  to  be  applied 
time  to  print  output 
2*PI 

rotor  forward  speed  (advance  velocity) 

resultant  velocity  at  blade  section  from  VU  and  VTHETA 

rotor  tip  speed  due  to  OMEGA 

velocity  comonent  at  blade  section  in  plane  of  rotation 

velocity  at  blade  section  normal  to  plane  of  rotation 

average  downwash  corresponding  to  TAVG 

initial  value  of  W  at  start  of  prediction  of  future  flapping 

dimensional  radius 

value  of  X  at  hub 
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Design  of  a  Digital  Controller 

As  with  the  numerical  approach,  a  dynamic  simulation  and  the  appropriate  feedback 
calculations  comprise  the  control  algorithm.  The  dynamic  simulation  predicts  the  rotor 
state  at  a  future  time  (usually  three  blade  revolutions  from  the  present),  and  the 
feedback  controller  determines  if  the  flapping  motion  will  be  excessive;  if  so,  the 
controller  then  automatically  provides  the  necessary  cyclic  control  step-input  to  limit  the 
flapping.  The  helicopter  dynamic  model  is  the  same  as  that  used  in  the  numerical 
control  approach. 

Implemented  as  a  FORTRAN  subroutine,  the  simulator  serves  the  dual  purpose  of 
generating  the  actual  rotor  motion,  and  of  predicting  the  future  motion  of  the  rotor, 
given  the  current  rotor  state  and  a  control  input.  In  the  actual  physical  system,  the 
rotor  state  and  control  input  would,  of  course,  be  determined  by  appropriate  sensors.  It 
is  envisioned  that  the  complete  controller  (including  the  dynamic  simulator)  would  be 
implemented  in  a  small  system  of  microprocessors. 

This  portion  of  the  report  describes  the  digital  controller,  presents  the  results  of 
several  simulations  to  test  its  performance,  and  concludes  with  a  feasibility  analysis  of 
implementation  using  dedicated  microprocessors. 

Construction  of  the  Closed-Loop  Control  System 

The  digital  controller  was  designed  using  largely  standard  techniques,  although  the 
feedback  gain  matrix  K  was  made  time-varying  in  order  to  give  better  performance  of 
the  system.  As  with  the  numerical  controller  in  the  earlier  part  of  this  report,  a 
simulator  was  used  to  generate  the  helicopter  rotor  flapping  motion. 

A  block  diagram  of  the  complete  system  appears  in  Fig.  19.  The  difference  between 
the  simulator  output  and  a  reference  signal,  which  is  the  physical  constraint  on  safe 
flapping  angle  P,  is  taken  as  the  feedback  signal.  The  non-linear  element  (proportional 
gain  with  deadband)  is  used  to  obtain  better  stability  properties  of  the  system.  A 
tracking  filter  is  utilized  to  depress  feedback  signal  noise  and  at  the  same  time  retain 
the  ability  to  track  varying  inputs.  A  first-order  hold  is  used  for  simulator  input  to 
achieve  better  tracking  ability  to  varying  inputs  while  zero-order  holds  are  used 
elsewhere  for  simplicity,  following  standard  design  practice. 


1. 


Design  of  Control  Elements 
Design  of  the  first-order  hold 


The  first-order  hold  is  constructed  with  reference  to  Fig.  20.  For  a  single  input- 
single  output  (SISO)  system,  the  hold  is  modeled  as: 

Xt)  -  *.-!>  *  •  HJ.  v,  <  « 

lk-l  ‘k-2 


where: 

y  =  quantity  being  sampled 
t  -  time 
k  «  time  index 
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This  SISO  hold  is  then  extended  to  the  multiple  input-multiple  output  (MIMO)  case. 

2.  Design  of  the  Deadband 

The  deadband  size  is  determined  by  system  performance  as  well  as  tracking 
performance.  The  system  performance  here  means  the  stability  and  transient 
characteristics  of  the  system.  From  the  stability  point  it  is  desirable  to  have  a 
relatively  large  deadband  size.  Alternatively,  from  the  perspective  of  tracking 
performance,  the  deadband  size  should  be  as  small  as  possible,  requiring  a  compromise. 

As  a  practical  consideration,  because  there  is  computing  error  in  the  simulator,  it  is 
undesirable  to  have  the  deadband  size  too  small.  In  the  actual  implementation  in  a 
helicopter,  sensor  noise  would  result  in  the  same  consideration.  For  the  given  helicopter 
(AH-1J),  the  deadband  A  is  chosen  as  0.2  degree  (as  shown  in  Fig.  21).  For  simplicity, 
a  and  P  are  both  chosen  as  45°  (i.e.,  the  nonlinear  element  has  unity  gain);  overall 
feedback  gain  is  adjusted  via  the  feedback  gain  matrix  K. 

3.  Design  of  the  feedback  gain  matrix. 

As  is  common  to  all  control  systems,  the  system  stability  is  critically  dependent  on 
the  value  of  feedback  gain.  Because  the  system  performance  is  not  easy  to  evaluate 
analytically  and  rotor  response  to  a  step  input  generally  takes  less  than  3  revolutions  to 
reach  its  steady  state,  it  is  sufficient  to  consider  a  first-order  simulator  for  the  present. 


As  expected,  an  improper  choice  of  feedback  gain  can  lead  the  system  to  diverge 
rapidly.  To  obtain  minimum  settling  time  of  the  system,  the  simulator  has  to  have  some 
overshoot  to  a  step  input.  This  requires  that  the  loop  gain  cannot  be  too  small.  On  the 
other  hand,  physical  considerations  do  not  allow  overshoot  in  the  system  (e.g.  an 
airframe  strike  could  be  the  result);  the  feedback  gain  must  not  be  so  large  that  the 
rotor  displays  an  oscillatory  transient.  To  resolve  this  contradiction,  a  varying  gain 
technique  is  used  so  that  the  gain  decreases  with  time  and  reaches  a  steady  state  value. 
The  initial  value  of  the  feedback  gain  can  be  chosen  so  that  the  simulator  prevents  an 
oscillatory  transient.  After  several  sampling  points,  the  gain  decreases  so  that  the 
transient  of  the  simulator  is  no  longer  oscillatory. 

(1)  Selection  of  initial  value  of  feedback  gain. 

Considering  the  closed  form  solution  of  rotor  flapping,  it  is  helpful  to  understand 
the  rotor  behavior.  From  the  closed  form  solutions,  we  have: 

P  -  -  Aj  cos  rf>  -  Bj  sini^ 

#0  Wf  3po  rrf4 

3e  “  A  ’  ~c6j  *  — S~ 

9Aj  rF  (A^  f4  -A  M  fj)  +  A  M  K.p 

mi - s - 


dA1 
sr 


1  dAi 

39a 


^  IV  (^12^4  “  +  ^14 

- S - 


aAl  K/,rr  -  1 

30^  “ - S - 


'2 
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®l  (^12^4"  AMf2)  +  ^14  ^-p  +  ®11^F^4 

m; - a 

where: 

0O  *  collective  pitch  angle 
61  *  lateral  cyclic  pitch  angle 
0j  ■  longitudinal  cyclic  pitch  angle 


A  -  (1  +  K/am)  (1  -  K/2rF)  +  KBf4rF  (Bu  +  k/au) 

Aj  «  longitudinal  flapping  angle 


Bj  -  lateral  flapping  angle 


Po  »  coning  angle 


3B 


-Pm  20 


B  -  tip  loss  factor 


r,-¥ 


f2  -  ^  +  Ba) 


,  B3  .  J*3 


fs  -  Ba  <-f-  ♦  -f-) 


*4*3 


fi  B* 


.  .  ,A*B3  m3 

^11  “  4  ("2”  '  T' 


.  8MB 

3(B3  - 
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b2  +  4 

o2  V 


C0  +  cT 


p  .  5.7  (%)4  .  j  /  MIF 


11  *  2  a2 

3(B2  - 

MIF  =  blade  moment  of  inertia  about  flapping  axis 

Kb  *  pitch-flap  coupling  coefficient. 

For  the  example  helicopter  (AH-1J), 

4  “  0'°  ’  ^  “  014 

3A,  3A, 


Obviously, 

dS 

wr  *  -  wr  sinV»  *  +  sin$ 

<Wj  Wj 

cosy  *  -  cosy 

Referring  to  Figs.  (22a-22b),  it  is  seen  that  the  closed  form  expressions  for 

apo  apo  3Aj  9Al  apx 

35^  ’  36£  ’  3§£  ’  ^  ^ 

are  close  to  the  values  calculated  numerically.  Similar  agreement  is  obtained  for  other 
values  of  the  parameters  o,  0O,  and  i>.  A  feedback  gain  of  1.0  gives  a  critically 
damped  system;  therefore,  the  initial  value  of  feedback  gain  should  be  chosen  larger  than 
1.0.  The  time-varying  gain  is  computed  as  follows: 

A  Norm  R  defined  for  a  vector  (n  x  1)  is: 


Norm  (U) 


t  e‘e  0 

L  0  e(  1)e  J 


This  measure  puts  some  weighting  on  each  element  of  U.  Note  that  Norm  (U)  3 
Euclidean  norm  of  U. 
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Let  U  be  a  10  x  1  vector  matrix,  and  define 


U(t,10)  -  E(t) 

“|U(t,l)|1/J  “ 

U(t,9)  =  E(t-1)  , 

U(t)  - 

|U(t,2)|1/2 

• 

• 

U(t,l)  =  E(t-9) 

».  |U(t,10)| 1/2  _ 

where  E  is  the  difference  signal  between  the  simulator  output  and  reference  signal,  as 
shown  in  Fig.  19. 

Then,  the  gain  is  deduced  using 


Gain  (t) 


Gain  (t-1) 


Norm  fU(t)] 
Norm  [U(t-l)J 


In  order  to  keep  the  tracking  ability,  a  lower  limit  is  set  on  gain.  In  the  case  studied, 
c  has  been  chosen  as  0.1. 


(2)  Selection  of  lower  limit  on  feedback  gain. 

The  choice  of  lower  limit  on  feedback  gain  is  made  on  the  basis  that  the  simulator 
should  not  be  oscillatory  after  some  time  t  and  the  settling  time  of  the  simulator  should 
be  small.  Obviously,  a  choice  of  dominant  damping  ratio  of  f  «  .707  (i.e.,  V2/2)  of 
the  simulator  dynamics  is  reinforced.  In  the  given  case,  the  lower  limit  is  chosen  as  0.6. 

4.  Design  of  the  Tracking  Filter 

For  the  SISO  system,  the  tracking  filter  is  shown  in  Fig  23(a),  and  the 
corresponding  root  locus  appears  in  Fig  23(b).  This  is  a  sampled-data  system,  and  the 
design  is  performed  accordingly  using  z-transforms.  Two  open-loop  poles  are  placed  at  z 
=  1,  so  that  the  filer  has  zero  steady  state  error  in  tracking  a  ramp  input. 

For  simplicity,  a  fixed  relation  between  a  and  k  is  chosen  as 


where 


k 


m 


2 


k  *  tracking  filter  gain 
a  -  open-loop  zero 

This  relation  has  been  determined  to  give  good  filtering  performance. 
The  filter  is  then; 


F(z)  -  2(a-l)  z  (z-a) 

v  '  (a-2)  z2  -  2a(a-2)  z-a 

To  illustrate  the  filtering  effect  of  F(z),  some  random  variations  (noise) 
superimposed  onto  a  deterministic  ramp  with  a  slope  of  0.S  are  generated  and  input  to 
the  system.  Referring  to  Figs.  (24a-24c),  the  tracking  results  are  given  for  different 
choices  of  a. 


It  is  seen  that  when  a  -  0.2,  the  filter  response  is  fast  enough  to  track  the  noise. 
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whereas  some  noise  depression  effect  is  achieved  when  a  ■  0.8.  As  the  filter  response 
speed  is  further  reduced,  the  filter  will  show  less  effect  to  noise  and  be  more  sluggish  in 
tracking  a  varying  true  signal. 

As  a  compromise,  a  =0.7  is  chosen.  This  SISO  filter  is  then  extended  to  our  MIMO 
system. 


Digital  Simulation  Results 

Simulations  were  done  with  the  sampling  period  being  the  time  needed  for  the  rotor 
blade  to  travel  30°. 

It  can  be  seen  from  Figs,  (25-26)  that  the  system  can  follow  a  step  input  reasonably 
with  no  overshoot  across  the  physical  Happing  angle  constraint  line  (the  dot-dash  line  in 
the  figures). 

From  Figs.  (27-28),  it  is  seen  that  there  is  a  steady-state  error  to  ramp  inputs  (i.e. 
the  flapping  actually  reaches  the  physical  limit).  This  steady-state  error  cannot  be 
eliminated  because  the  control  is  based  on  an  imperfect  prediction  of  future  system 
response.  Thus,  no  integral  component  in  the  loop  can  be  used  to  eliminate  the  steady 
state  error.  Fortunately,  this  error  is  not  very  big  and  can  be  reduced  if  a  higher-order 
hold  is  used  for  control  inputs.  Practically,  a  physical  constraint  on  safe  flapping  angle 
can  be  set  so  that  the  ramp-tracking  error  of  the  closed-loop  system  may  be  allowed 
with  no  danger. 


Implementation  Requirements 

Figure  29  is  a  block  diagram  of  a  typical  microprocessor  arrangement  used  in  a 
control  system.  The  input  device  includes  sensors  and  an  A/D  converter,  while  the 
output  device  includes  a  D/A  converter.  The  filter  (digital  to  analog)  is  usually  employed 
to  reject  noise  from  the  senors  and  to  smooth  the  control  output.  Coprocessors  may  be 
utilized  to  reduce  the  complexity  of  the  software  and/or  increase  the  processing  speed 
of  the  system  through  parallelism.  For  the  helicopter  blade  flapping  problem,  the  input- 
output  variables  are  shown  in  Table  2. 

1.  Sampling  Rates 

For  digital  control,  sampling  is  necessary  and  the  choice  of  sampling  rate  is  crucial  to 
the  performance  of  the  system.  The  sampling  of  the  measurement  signals  ft  and 

ft  must  be  based  upon  azimuthal  angle  rather  than  being  time-based,  because 

in  this  application,  ft  and  ft  are  periodic  with  respect  to  i>.  It  was  determined 
that  a  desirable  sampling  rate  would  be  at  least  5  times  per  revolution.  Clearly,  the 
value  of  the  sampling  interval  must  be  less  than  72  degrees.  For  a  main-rotor  tip- 
speed  of  approximately  740  ft/sec,  and  a  blade  diameter  of  approximately  44  ft,  the 
sampling  period  (in  time)  must  by  less  than  75  ms.  Such  a  rate  is  readily  achieved  with 
existing  technology. 
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2.  Storage  Requirements 

There  are  approximately  150  variables  to  be  stored  in  24-bit  RAM  for  the 
simulation.  Information  on  CL  and  CD  and  all  other  helicopter  parameters  would  be 
stored  in  approximately  100  words  of  24-bit  ROM.  In  addition,  the  trigonometric 
function  "sin”  would  be  stored  in  half-degree  increments  from  0  to  90  degrees  in  a  look¬ 
up  table  (for  3-digit  precision,  this  would  require  180  words  of  16-bit  ROM).  The 
controller  requires  another  50  variables  in  24-bit  RAM.  The  computational  processor 
would  need  another  50  words  of  24-bit  RAM.  Therefore,  the  total  memory  requirement 
for  implementation  is 

RAM:  750  bytes 
ROM:  660  bytes 

which  is  a  relatively  trivial  amount  of  memory  for  current  technology.  For  example,  the 
AMD  80C51  CMOS  single-chip  controller  has  4K  bytes  of  on-chip  ROM  and  128  bytes  of 
on-chip  RAM;  therefore,  only  IK  bytes  of  external  RAM  are  needed  for  this  application. 

3.  Microprocessor  Speed  Requirement 

To  calculate  the  arithmetic  operations  done  in  one  sampling  period,  the  software 
flow-chart  of  the  system  to  be  implemented  is  shown  in  Fig.  29.  The  operations  in  each 
stage  are  arranged  in  Table  3.  If  parallel  processing  is  employed  for  the  inner  loop  in 
Fig.  29,  the  necessary  operations  are  as  shown  in  Table  4.  To  execute  all  of  these 
instructions,  an  AMD  80C51  would  require  about  900  ms  with  its  CPU  operating  at  12 
MHz  (working  with  24-bit  data). 
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Conclusions  and  Recommendations 

It  appears  to  be  feasible  to  limit  the  amount  of  flapping  of  a  helicopter  rotor  blade 
by  means  of  a  numerical  controller  which  continuously  predicts  the  future  flapping  of  the 
rotor  given  the  instantaneous  state  of  the  rotor  and  control  input.  Admittedly,  this 
present  study  is  limited  in  scope  but  enough  has  been  learned  about  the  response  of  the 
rotor  to  the  proposed  feedback  system  to  justify  proceeding  further. 

One  of  the  major  unanswered  questions  is  that  of  the  speed  with  which  a 
microprocessor  might  do  the  calculations  of  the  future  rotor  flapping.  As  an  order  of 
magnitude,  it  appears  as  if  the  microprocessor  might  be  called  upon  to  calculate  the 
rotor  flapping  two  revolutions  ahead  in  the  time  it  takes  the  rotor  to  move  approximately 
0.2  of  a  revolution,  or  around  40  milliseconds.  This  is  certainly  not  out  of  the  question, 
particularly  considering  the  possibility  that  one  might  be  able  to  use  closed-form 
equations  to  predict  the  instantaneous  aerodynamic  moment  on  the  blade.  In  this  case, 
numerical  integrations  need  only  be  done  with  respect  to  the  azimuth  angle. 

For  the  digital  controller,  the  time  required  for  the  microprocessor  appears  at  first 
to  be  prohibitively  large,  but  the  implementation  estimates  were  made  using  a  relatively 
inexpensive  "off  the  shelf"  microprocessor.  A  custom-made  chip  could  easily  complete 
the  required  calculations  in  one-tenth  to  one-twentieth  of  the  time  estimated  here.  At 
the  current  rate  of  increase  in  processor  speeds,  it  appears  likely  that  a  standard 
counterpart  to  the  microprocessor  analyzed  here  (AMD  80C51),  with  ten  to  twenty  times 
its  speed,  will  be  available  within  a  year  or  two.  As  with  the  numerical  controller, 
further  decreases  in  required  computations  might  be  achieved  using  closed-form  equations 
for  the  blade-aerodynamics. 

Unlike  a  stability  augmentation  system,  the  authority  of  the  system  presented  here 
will  probably  have  to  be  relatively  high.  However,  this  system  will  only  provide  a  signal 
to  cyclic  control  actuators  at  those  rare  times  when  the  rotor  is  disturbed  to  such  a 
degree  as  to  result  in  excessive  flapping. 

If  further  work  is  pursued  on  this  system,  the  following  recommendations  are  made. 

a.  The  equations  of  motion  of  the  airframe  should  be  included. 

b.  Considerably  more  cases  should  be  examined,  including  trim  at  lower  airspeeds 
with  combinations  of  collective  and  cyclic  control  inputs. 

c.  The  response  of  teetoring  rotors  should  be  compared  with  articulated  rotors. 

d.  The  time  required  for  a  microprocessor  to  do  the  required  calculations  should 
be  studied  realizing  that,  in  the  operating  system,  the  controller  will  only  have 

to  predict  the  future  flapping  and  not  the  continuous  flapping  as  is  done  here 
for  the  simulation. 
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Nomenclature  for  Digital  Control  Program  (DCFLAP) 

Note:  Because  the  digital  control  program  uses  the  same  helicopter  simulation  code 
as  the  numerical  control  program,  many  of  the  variables  are  common  to  both.  This  table 
generally  lists  those  variables  that  appear  only  in  the  digital  program  (DCFLAP).  Refer 
to  the  Nomenclature  for  the  Numerical  Control  Program  for  other  variable  definitions. 

Variable  Name  Definition 

All  temporary  dummy  variable 

A12  temporary  dummy  variable 

A13  temporary  dummy  variable 

A14  temporary  dummy  variable 

AIRFOIL  subroutine  to  calculate  section  lift  and  drag  coefficient  CL  &  CD 
ALPHA  initial  trim  angle  of  attack  of  disk  plane 

ALPHAD  trim  angle  of  attack  of  disk  plane  in  degrees 

ANORM  subroutine  to  calculate  ANORM(U) 

ANORM(U)  a  special  norm  of  a  vector  with  data  weighting 
Bll  temporary  dummy  variable 

BEEP  logical  variable,  signaling  occurrence  of  excessive  flapping 

BET  value  of  BETA  in  degrees  for  printout  or  writing  to  initial  condition 

table 

BET  10  initial  trim  flapping  angle  in  radians  at  PSIBD 

BET10D  temporary  dummy  variable 

BETID  1  temporary  dummy  variable  in  interpolation  of  initial  condition  table 

BET1D2  temporary  dummy  variable  in  interpolation  of  initial  condition  table 

BET20  first  derivative  of  flapping  angle  corresponding  to  BET10 

BET20D  temporary  dummy  variable 

BET2D1  temporary  dummy  variable  in  interpolation  of  initial  condition  table 

BET2D2  temporary  dummy  variable  in  interpolation  of  initial  condition  table 

value  of  BET  AD  in  degs/sec  for  printout  or  writing  to  initial  condition 
table 


BETD 


BWFLAP 

CDELAY 

CNTRL 

CNTRL1 

CNTRL2 

CRATEA 

CRATEO 

DO 

D1 

D2 

DALPHA 

DAMP 

DELAY 

DENOM 

DETECT 

El 

E2 

FI 

F2 

F3 

F4 

FDBACK 

GAIN 

IK  TEN 
KTEN 
NBET 
NENTRY 


subroutine  to  simulate  the  rotor  blade  flapping  motion 
temporary  dummy  variable 
subroutine  formulating  control  law 

feedback  signal  of  lateral  cyclic  pitch  in  degrees  for  printout 

feedback  signal  of  longitudinal  cyclic  pitch  in  degrees  for  printout 

rate  of  increase  of  angle  of  attack  of  disk  plane 

rate  of  increase  of  collective  pitch 

perturbation  in  collective  pitch 

perturbation  in  lateral  cyclic  pitch 

perturbation  in  logitudinal  cyclic  pitch 

perturbation  in  angle  of  attack  of  disk  plane 

damping  ratio  of  rotor  blade  flapping 

phase  angle  delay  of  rotor  blade  flapping 

temporary  dummy  variable 

logic  switch,  initiating  perturbation  if  NREV  reaches  PERTB 

THETA  1  in  degrees  for  printout 

THETA2  in  degrees  for  printout 

temporary  dummy  variable 

temporary  dummy  variable 

temporary  dummy  variable 

temporary  dummy  variable 

subroutine  to  determine  the  amount  of  feedback 

magnitude  of  feedback  gain  vector,  see  Fig.  19  (the  system  block 
diagram) 

logic  switch  for  sampling 

sampling  is  done  once  every  KTEN  integration  points 
number  of  points  per  revolution  in  initial  condition  table 
entry  point  in  initial  condition  table 
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NREV 

OVER 

OVERO 

PSIOO 

PSI1ST 

PSIBD 

PSIBD1 

PSID 

REPLY 

SOS 

STEP 

T1 

T2 

T3 

T4 

THETO 

THETOD 

THET1 

THET1D 

THET2 

THET2D 

THRST1 

THRST? 

TRIM 

TTAVG 

U0(I) 


azimuth  angle  PSI  in  terms  of  number  of  revolutions  from  the  start  of 
simulation 

amount  of  excessive  flapping 
maximum  amount  of  excessive  flapping 
temporary  dummy  variable 

temporary  dummy  variable  in  preparing  the  initial  condition  table 

azimuth  angle  in  degrees  when  simulation  begins 

temporary  dummy  variable  in  interpolation  of  initial  condition  table 

value  of  PSI  in  degrees  for  printout  or  writing  to  initial  condition 
table 

logic  switch,  user  keyboard  input 

azimuth  position  at  which  maximum  excessive  flapping  occurs 

step  size  of  PSI  in  degrees  in  the  initial  condition  table 

temporary  dummy  variable 

temporary  dummy  variable 

temporary  dummy  variable 

temporary  dummy  variable 

initial  trim  collective  pitch 

trim  collective  pitch  in  degrees 

initial  trim  lateral  pitch 

trim  lateral  cyclic  pitch  in  degrees 

initial  trim  longitudinal  pitch 

trim  longitudinal  cyclic  pitch  in  degrees 

temporary  dummy  variable 

temporary  dummy  variable 

logical  variable,specifying  the  status  of  existence  of  initial  condition  table 

temporary  dummy  variable 

similar  to  U(I)  but  used  as  a  temporary  one 


UTHETl(I)  temporary  dummy  array,  used  in  forming  first-order  hold 
UTHET2(I)  temporary  dummy  array,  used  in  forming  first-order  hold 
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U(I)  vector  windowing  the  flapping  angle 

VAR1  temporary  dummy  variable 

VAR2  temporary  dummy  variable 

VAR3  temporary  dummy  variable 

VNUL  temporary  dummy  variable 

WNREL  relative  natural  frequency  of  rotor  blade  flapping  w.r.t  OMEGA 
XA1  temporary  dummy  variable 

XB1  temporary  dummy  variable 


31 


TABLE  1 


Parameters  Corresponding  to  Figures  13  through  17 

CASE 

DELFB1 

DELFB2 

FB1LTM  FB2LTM 

NPRED 

FPRED 

1* 

2 

2 

5 

5 

1 

.2 

2* 

2 

2 

5 

5 

1 

.2 

3 

2 

2 

5 

5 

1 

.2 

4 

2 

2 

5 

5 

2 

.2 

5 

2 

2 

5 

5 

3 

.2 

6 

4 

4 

8 

8 

3 

.1 

7 

4 

4 

8 

8 

3 

.2 

8 

4 

4 

8 

8 

3 

.4 

9 

4 

4 

8 

8 

2 

.2 

10 

3 

3 

8 

8 

2 

.2 

11 

2 

2 

8 

8 

2 

.2 

12 

2 

2 

6 

6 

2 

.2 

13 

2 

2 

7 

7 

2 

.2 

14 

4 

4 

8 

8 

3 

.1 

15** 

2 

2 

8 

8 

2 

.2 

16 

4 

4 

7 

7 

1 

.2 

17 

4 

4 

7 

7 

2 

.2 

18 

4 

4 

7 

7 

3 

.2 

FOR  ALL  CASES 

BETLIM 

*  8 

(•BETLIM  - 

50) 

TOFF 

-  10 

PRT 

-  3 

NMAX 

«  10 

CRATE  1 

=  100 

(**CRATE1 

-  200) 

CRATE2 

=  100 

(**CRATE2 

=  200) 
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TABLE  2 

Input-Output  Variables  for  Digital  Controller 


VARIABLE 

INPUT 

OUTPUT 

Downwash,  w 

X 

Thrust,  T 

X 

Azimuth  position,  ^ 

X 

Flapping  angle,  0 

X 

• 

Flapping  rate,  0 

X 

Forward  speed,  v 

X 

Angle  of  attack,  « 

X 

Collective,  60 

X 

Lateral  cyclic, 

X 

X 

Longitudinal  cyclic,  02 

X 

X 

Air  density,  p 

X 

4 
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TABLE  3 

Operations  in  Each  Stage 


STAGE 

ADD 

SUB 

MIZL 

DIV 

REF. 

1 

4 

6 

29 

12 

9 

2 

2 

3 

11 

5 

9 

3 

15 

15 

35 

12 

17 

4 

69 

55 

105 

22 

34 

STAGE 

CMP 

BRANCH  ABS 

SORT 

OPS 

1 

0 

0 

0 

1 

0 

2 

0 

0 

0 

0 

0 

3 

111 

3 

4 

1 

4 

4 

10 

7 

5 

0 

3 
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add 

SUB 

MUL 

rnv 

MEM. 

EEEi 

90,000 

7,800 

18,000 

5,100 

6,900 

CMP 

BRANCH  ABS 

S-QRI 

BINARY 

OPS 

12,100 

1,000 

900 

200 

700 

FIGURE  1 

LEFT  SIDE  VIEW  OF  A  ROTOR  SHOWING  THE 
ANGLE  OF  ATTACK  AND  FREESTREAM  VELOCITY  COMPONENTS 
PARALLEL  AND  NORMAL  TO  THE  DISC  PLANE 


FIGURE  2 

TOP  VIEW  OF  THE  DISC  PLANE  SHOWING  VELOCITY 
COMPONENTS  PARALLEL  AND  NORMAL  TO  THE  BLADE 
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view  in  plane  of 
blade  and  shaft  axis 


FIGURE  3 

BLADE  -  SHAFT  AXIS  PLANE  SHOWING  VELOCITY  COMPONENTS 
RESULTING  FROM  FLAPPING  AND  DOWNWASH 


FIGURE  4 

BLADE  AIRFOIL  SECTION  SHOWING  THE  TANGENTIAL 
AND  NORMAL  VELOCITIES  TO  THE  DISC  PLANE  WHICH 
DETERMINE  THE  SECTION  ANGLE  OF  ATTACK 


I 

o 

CO 


I 

o 

CM 


r  r 

©  o 


T 


saaaoaa  ‘ONidd\ru  navne  «oloh 


Ot- 


39 


•21  M 


00 


angle,  thousands 
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NOTE:  FIGURE  18  INTENTIONALLY  OMITTED. 


Figure  19  Block  diagram  of  the  digital  control  system. 


g  angle  p 
Odeg. 


i  angle  p 
0  deg. 
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Figure  23a.  Block  diagram  of  the  tracking  filter. 


Figure  24a  Filter  output  for 


a 

°  3 
a  A 


w  9 

a-f 
a  § 
■a  s 

ii 

I  o 


(b^jt3t\  ;nd;no  J»ma 


T 


T 


T 


T 

o 


r 

o  o 

cd  « 


o  o  o 

W  CM  *H 


(safari  ^nd^o  JtaTTU 


Figure  24c  Filter  output  for  a  -  0.99 
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Figure  26c  Longitudinal  cyclic  pitch  02  for  a  s 
of  A02  =  -6  deg.  at  NREV  =  2 

(with  digital  controller) 


2.4 
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NTV9 


Figure  26d  Feedback  gain  for  a  step  input 
of  A02  =  -6  deg.  at  NREV  =  2 
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(  gsp)  IQ 


0,-5  deg/sec  at  NRBV  - 

(with  di|iul  controller) 
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bnvf) 


Figure  27d  Feedback  gain  for  ramp  input  of 
§4-5  deg/sec  at  NRBV  -  2 


Ua  coo 


Figure  28b  Lateral  cyclic  pitch  0 1  for  a  ramp  input 
8j  -  5  deg/sec  at  NREV  -  2 

(with  digital  controller) 
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(•8sp)  *0 


Figure  28c  Longitudinal  cyclic  pitch  02  for  a  ramp  input  of 
02  =  5  deg/sec  at  NREV  =  2 

(with  digital  controller) 


2.6 


hcrvo 


Figure  28d  Feedback  gain  for  a  ramp  input  of 
CL  »  5  deg/sec  at  NRBV  «  2 
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Figure  29  Block  diagram  of  the  microprocessor  controller. 


30  Software  flowchart  for  analyzing  microprocessor  operations. 
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FORTRAN  listing  for  PROGRAM  FLAP 
(numerical  controller) 
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PROGRAM  FLAP 

COMMON/ COM1/THE1LM , THE2LM , KNTRL , TCON , TOFF , THE1I , THE2I , SWITCH 
COMMON/ COM2/DELFB1 , DELFB2 , FB , PI , CRATE1 , CRATE2 , FDBK1 , FDBK2 
COMMON/COM3/CLI (20) ,ALPHAL(20) ,CDI(33) ,ALPHAD(33) , ALPHAB, CL, CD 
C0MM0N/C0M4/0MEGA , V , ALPHA , W , BETAD , RHO , DELPS I , DELR , T , MI F , MW , D 
COMMON/ COM5/PS I , EPS , CO , CT , THETAO , THETAT , THETAl , THETA2 , KBETA , BETA 
C0MM0N/C0M6 /FB1LIM , FB2LIM 
REAL  MIF, MW, KBETA, N,NPRED,NMAX 
OPEN (UNIT-1 , FILE- ' CL . DAT ' ) 

DO  100  1-1,20 

READ(1 ,*)  ALPHAL(I) ,CLI(I) 

100  CONTINUE 

CLOSE  (UNIT-1) 

OPEN (UNIT-2 , FILE- ' CD . DAT ' ) 

DO  200  1-1,33 

READ(2 , *)  ALPHAD(I) ,CDI(I) 

200  CONTINUE 

CLOSE (UNIT-2) 

PI-3 . 14159 
TWOPI-2 .*PI 
DTR— PI/180 . 

DO  300  1-1,20 
ALPHAL ( I ) -ALPHAL ( I ) *DTR 
300  CONTINUE 

DO  400  1-1,33 
ALPHAD ( I ) -ALPHAD ( I ) *DTR 
400  CONTINUE 

WRITE(* ,*) 'THIS  PROGRAM  PREDICTS  THE  FLAPPING  MOTION  FOR  A' 

WRITE (*,*) 'HELICOPTER  ROTOR  BLADE  IN  FORWARD  FLIGHT  WITH  PILOT' 
WRITE (*,*)' INPUT  AND  FEEDBACK  CONTROL.' 

WRITE(*,*) 

WRITE (*,*) 

WRITE (*,*) 'INPUT  IDENTIFYING  CASE  NUMBER' 

READ(*,*)CASE 

WRITE(*,*)' INPUT  THE  MAXIMUM  FLAPPING  ANGLE  IN  DEGREES  WHICH  ' 
WRITE (*,*) 'WILL  BE  ALLOWED' 

READ(*,*)BETLIM 

WRITE(*,*) 'INPUT  THE  INCREMENTAL  CORRECTION  TO  LATERAL  CYCLIC' 
WRITE(*,*)' PITCH  CONTROL  WHICH  IS  ADDED  EACH  TIME  BETA  IS' 
WRITE(*,*) 'PREDICTED  TO  EXCEED  THE  LIMIT.  IT  IS  ALSO  SUBTRACTED' 
WRITE(* ,*) ' EACH  TIME  BETA  IS  PREDICTED  NOT  TO  EXCEED  THE  LIMIT’ 
WRITE(*,*) ' IF  ANY  FEEDBACK  IS  BEING  APPLIED.' 

READ(*,*)DELFB1 

WRITE (*,*) 'SIMILARLY,  INPUT  THE  INCREMENTAL  CORRECTION  TO' 
WRITE(*,*) 'THE  LONGITUDINAL  CYCLIC' 

READ(*,*)DELFB2 

WRITE(*,*)' INPUT  THE  LIMIT  ON  FEEDBACK  TO  LATERAL  CYCLIC , DEGS . ' 
READ ( * , * ) FB1LIM 

WRITE(*,*)' INPUT  THE  LIMIT  ON  FEEDBACK  TO  LONG.  CYCLIC,  DEGS.' 
R£AD(*,*)FB2LIM 

WRITE(*,*) 'INPUT  NUMBER  OF  ROTATIONS  BEFORE  INITIATING  CONTROL' 
READ(*,*)N 

WRITE (*,*) 'INPUT  NUMBER  OF  ROTATIONS  TO  PREDICT  AHEAD' 
READ(*,*)NPRED 
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PROGRAM  FLAP 

COMMON/COM1/THE1LM , THE2LM , KNTRL , TCON , TOFF , THE1I , THE2I , SWITCH 
COMMON/COM2/DELFB1 , DELFB2 , FB , PI , CRATE1 , CRATE2 , FDBK1 , FDBK2 
COMMON/COM3/CLI (20) ,ALPHAL(20) ,CDI(33) ,ALPHAD(33) , ALPHAB , CL , CD 
COMMON/ C0M4/0MEGA , V , ALPHA , W , BETAD , RHO , DELPSI , DELR , T , MIF , MW , D 
COMMON/COM5/PSI , EPS , CO , CT , THETAO , THETAT , THETAl , THETA2 , KBETA , BETA 
COMMON/COM6 /FB1LIM , FB2LIM 
REAL  MIF , MW , KBETA , N , NPRED , NMAX 
OPEN (UNIT-1 , FILE- ' CL . DAT ' ) 

DO  100  1-1,20 

READ(1 ,*)  ALPHAL(I) , CLI(I) 

100  CONTINUE 

CLOSE  (UNIT-1) 

OPEN (UNIT-2 , FI LE- ' CD . DAT ' ) 

DO  200  1-1,33 

READ(2 ,*)  ALPHAD (I) ,CDI(I) 

200  CONTINUE 

CLOSE (UNIT-2) 

PI-3.14159 
TWOPI— 2.*PI 
DTR— PI/180. 

DO  300  1-1,20 
ALPHAL ( I ) -ALPHAL ( I ) *DTR 
300  CONTINUE 

DO  400  1-1,33 
ALPHAD ( I ) -ALPHAD ( I ) *DTR 
400  CONTINUE 

WRITE (*,*) 'THIS  PROGRAM  PREDICTS  THE  FLAPPING  MOTION  FOR  A' 

WRITE (*,*) 'HELICOPTER  ROTOR  BLADE  IN  FORWARD  FLIGHT  WITH  PILOT' 
WRITE (*,*)' INPUT  AND  FEEDBACK  CONTROL.' 

WRITE(*, *) 

WRITE(*,*) 

WRITE (*,*) 'INPUT  IDENTIFYING  CASE  NUMBER' 

READ (*,*) CASE 

WRITE (*,*) 'INPUT  THE  MAXIMUM  FLAPPING  ANGLE  IN  DEGREES  WHICH  ' 
WRITE (*,*) 'WILL  BE  ALLOWED' 

READ(*,*)BETLIM 

WRITE(*,*)' INPUT  THE  INCREMENTAL  CORRECTION  TO  LATERAL  CYCLIC' 
WRITE (*,*) 'PITCH  CONTROL  WHICH  IS  ADDED  EACH  TIME  BETA  IS' 
WRITE(*,*) 'PREDICTED  TO  EXCEED  THE  LIMIT.  IT  IS  ALSO  SUBTRACTED' 
WRITE(*,*) 'EACH  TIME  BETA  IS  PREDICTED  NOT  TO  EXCEED  THE  LIMIT' 
WRITE(*,*) ' IF  ANY  FEEDBACK  IS  BEING  APPLIED.' 

READ(*,*)DELFB1 

WRITE(*,*) 'SIMILARLY,  INPUT  THE  INCREMENTAL  CORRECTION  TO' 
WRITE(*,*) 'THE  LONGITUDINAL  CYCLIC' 

READ(*,*)DELFB2 

WRITE(*,*)' INPUT  THE  LIMIT  ON  FEEDBACK  TO  LATERAL  CYCLIC  DEGS.' 
READ(*,*)FB1LIM 

WRITE(*,*)' INPUT  THE  LIMIT  ON  FEEDBACK  TO  LONG.  CYCLIC,  DEGS.' 
READ(*,*)FB2LIM 

WRITE (*,*) 'INPUT  NUMBER  OF  ROTATIONS  BEFORE  INITIATING  CONTROL' 
READ(*,*)N 

WRITE(*,*)' INPUT  NUMBER  OF  ROTATIONS  TO  PREDICT  AHEAD’ 

READ (*,*) NPRED 
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WRITE (*,*) 'INPUT  FRACTION  OF  REV  REQUIRED  TO  PREDICT' 

WRITE (*,*) 'BLADE  MOTION  NPRED  REVOLUTIONS  AHEAD' 

READ(*,*)FPRED 

WRITE (*,*) 'ONLY  CYCLIC  CONTROLS  WILL  BE  APPLIED.  EACH  CONTROL' 
WRITE (*,*) 'WILL  BE  INCREASED  LINEARLY  UP  TO  A  MAX  VALUE.' 

WRITE (*,*) 'INPUT:  RATE  OF  INCREASE  OF  LATERAL  CYCLIC,  DEGS /SEC' 
WRITE(*,*) '  RATE  OF  INCREASE  OF  LONGITUDINAL  CYC.DEGS/SEC' 

READ (*,*) CRATE 1 , CRATE2 

WRITE (*,*) 'INPUT:  MAX  INCREMENTAL  VALUE  FOR  LATERAL  CONTROL,  DEG' 
WRITE(*,*) '  MAX  INCREMENTAL  VALUE  FOR  LONG.  CONTROL,  DEGS' 

READ ( * , * ) THE1LM , THE2LM 

WRITE (*,*) 'INPUT:  LENGTH  OF  TIME  FOR  CONTROL  TO  BE  APPLIED, SECS' 
WRITE(*,*) '  INPUT  NUMBER  OF  CALCULATIONS  BETWEEN  PRINTOUTS' 

READ (*,*) TOFF, PRT 

WRITE (*,*) 'INPUT  MAX  NUMBER  OF  REVS  FOR  RUN' 

READ(*,*)NMAX 

OPEN (UNIT-5 , FILE- ' FLP5 . DAT ' ) 

READ  (5,*)D,C0,CT,EPS 
CLOSE (UNIT-5) 

C  D-ROTOR  DIAMETER, FT. 

C  CO-ROOT  CHORD  AT  R-0,  FT. 

C  CT-TIP  CHORD,  FT. 

C  EPS-DIMENSIONLESS  HINGE  OFFSET 

OPEN (UNIT-6 , FILE- ' FLP6 . DAT ' ) 

READ  (6 ,*)KBETA,THETAT,MW,MIF,V 
CLOSE (UNIT-6) 

C  KBETA-PITCH-FLAP  COUPLING 

C  THETAT-TOTAL  TWIST  FROM  ROOT  TO  TIP  IN  DEGS,  NEGATIVE  FOR  WASHOUT 
C  MW-BLADE  WEIGHT  MOMENT  ABOUT  FLAPPING  HINGE  IN  FOOT-POUNDS 
C  MIF-BLADE  MASS  MOMENT  OF  INERTIA  ABOUT  FLAPPING  AXIS 
C  V-FORWARD  VELOCITY  IN  KTS 

OPEN (UNIT— 7 , FILE- ' FLP7 . DAT ' ) 

READ  ( 7 , *) VT , THETAO , THETAl , THETA2 , ALPHA , Al , B1 , BETAQ 
CLOSE(UNIT-7) 

C  VT-TIP  SPEED  DUE  TO  ROTATION  IN  FPS 

C  THETAO-INITIAL  TRIM  COLLECTIVE  PITCH,  DEGS 

C  THETAl «INITIAL  TRIM  LATERAL  CYCLIC  PITCH,  DEGS. 

C  THETA2-INITIAL  TRIM  LONGITUDINAL  CYCLIC  PITCH,  DEGS. 

C  ALPHA-INITIAL  TRIM  DISC  PLANE  ANGLE  OF  ATTACK,  DEGS. 

C  Al-LONGITUDINAL  FLAPPING  IN  DEGREES 

C  BI-LATERAL  FLAPPING  IN  DEGREES 

C  BETAO-CONING  ANGLE  IN  DEGREES 

OPEN (UNIT-8 , FILE- ' FLP8 . DAT ' ) 

READ  (8,*)DELX,DELPSI,RHO,TAVG 
CLOSE (UNIT-8) 

WRITE(*, 509) 

WRITE ( * , 50  3 ) CASE , BETLIM 
WRITE ( * , 504 ) DELFB 1 , DELFB2 
WRITE(*,505)FB1LIM,FB2LIM 
WRITE (* , 506)CRATE1 , CRATE2 
WRITE ( * , 507 ) THE1LM , THE2LM 
WRITE(* , 508 )TOFF , N , NPRED 

503  FORMAT ( '  ' , 'CASE-' ,F9.3, '  BETLIM- ', F9 . 3 ) 

504  FORMAT ( '  ' , ' DELFB1-' , F9 . 3 , '  DELFB2-' , F9 . 3) 
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505 

506 

507 

508 

509 


501 

502 
C 

C 

C 

C 

C 

C 


C 


FORMAT (' 
FORMAT (' 
FORMAT (' 
FORMAT ( ' 
.'ORMAT  ( '  1 


,  ' FB1LIM-' , F9 . 3 , ' 
, ' CRATE1- ' , F9 . 3 , ' 
, 'THE1LM-' , F9 . 3 , ' 
,  'TOFF-' , F9 . 3 , ' 

) 


FB2LIM- ' , F9 . 3 ) 

CRATE2- ' , F9 . 3) 

THE2LM-' , F9.3) 

N- ' , F9 . 3 , '  NPRED— ' , F9 . 3) 


WRITE (*,*) 'PROGRAM  IS  HALTED  FOR  OPPORTUNITY  TO  PRINT  SCREEN' 
WRITE (*,*)' INPUT  1  (OR  ANYTHING)  TO  CONTINUE' 

READ (*,*) DUMMY 


OPEN (UNIT-20 , FILE- ' FLAP . DAT ' , STATUS- ' NEW' ) 

WRITE ( 20 , 501) CASE , BETLIM , DELFB1 , DELFB2 , FB1LIM , FB2LIM 
WRITE ( 20 , 502 ) N , NPRED , FPRED , CRATE1 , CRATE2 , THE1LM , THE2LM , TOFF 
FORMAT(6F9 . 3) 

FORMAT (8F9 . 3) 

DELX-INCREMENT  IN  DIMENSIONLESS  BLADE  RADIUS  FOR  NUMERICAL 
INTEGRATION  OF  THRUST  WITH  RADIUS 
DELPS I-INCREMENT  IN  AZIMUTH  ANGLE,  DEG3 ,  FOR  NUMERICAL 
INTEGRATION  OF  BLADE  MOTION  WITH  PSI 
RHO-AIR  MASS  DENSITY,  SLUGS/CU.  FT. 

TAVG-INITIAL  AVERAGE  THRUST  OF  ROTOR  IN  ONE  REVOLUTION 
AREA— PI*D*D/4 . 

V-V*l . 69 

BETLIM-BETLIM*DTR 

DELFB1-DELFB1*DTR 

DELFB2-DELFB2*DTR 

FB1LIM-FB1LIM*DTR 

FB2LIM-FB2LIM*DTR 

THETAO-THETAO*DTR 

THETAl-THETAl*DTR 

THETA2-THETA2*DTR 

THETAT-THETAT*DTR 

Al-Al*DTR 

B1-B1*DTR 

BETA0-BETA0*DTR 

THE1I-THETA1 

THE2I-THETA2 

CRATE1-CRATE1*DTR 

CRATE2— CRATE2*DTR 

THE1LM-THE1LM*DTR 

THE2LM-THE2LM*DTR 

ALPHA-ALPHA*DTR 

OMEGA-VT/D*2 

DELPS I-DELPS I*DTR 

DELT-DELPS I/OMEGA 

DELTPR-DELT*PRT 

PSIMAX-NMAX*TWOPI 


W-0.0 

DELR-DELX*D/2 
TPRED-FPRED*TWOPI/ .MEGA 
INITIALIZATION  BEGINS 
TPRINT-0 . 0 
TIME-0. 

PSI-0.0 
PSII-0.0 
BETA-BETAO - Al 
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BETAD— 0MEGA*B1 

SWTCH2-0 . 0 

FB-0.0 

DELT1-0 . 0 

KNTRL-1 

TCALOO .  0 

KNTRL-1 

PS ILIM-N*TWOPI 

TCONI-PS ILIM/OMEGA 

CALL  DWNWSH (V,TAVG,RHO, AREA, Al, ALPHA, W) 

1  CONTINUE 

IF(KNTRL. EQ. 1)G0  TO  2 
10  CONTINUE 

TCON— TIME-TCONI 

IF( PSI . GT . PSIMAX) GO  TO  5000 

CALL  SUBCNTRL 

2  CONTINUE 
CALL  SUBFLAP 

3  CONTINUE 
DELT2-T 

TCALC-TCALC+ (DELT1+DELT2 ) /2*DELPSI 
DELT1— DELT2 
TIME-TIME+DELT 
PSI— PSI+DELPSI 

4  CONTINUE 
IF(KNTRL.EQ.2)GO  TO  8 
IF(KNTRL.EQ.3)GO  TO  11 

5  CONTINUE 

IF(PSI.GE.PSII+TWOPI)GO  TO  19 
IF(PSI.GE.PSILIM)GO  TO  18 

6  CONTINUE 
IF(KNTRL.EQ. 2)GO  TO  7 
IF(KNTRL.EQ.3)GO  TO  7 

IF ( TIME . LT . TPRINT ) GO  TO  7 

TH1DG— THETAl/DTR 

TH2DG-THETA2/DTR 

REV— PSI/TWOPI 

BETADG-BETA/DTR 

FDBK1G-FDBK1/DTR 

FDBK2G— FDBK2/DTR 

WRITE(*,*) 'TIME-' .TIME, '  REV-* .kjlv ,  •  BETA-' .BETADG 
WRITE ( 20 , 500 ) TIME , TH1DG , TH2DG , FDBK1G , FDBK2G , REV , BETADG 
500  FORMAT(7F8.2) 

IF(PSI .GE. PSIMAX)GO  TO  5000 
TPRINT-TPRINT+DELTPR 

7  CONTINUE 
IF(KNTRL.EQ.l)GO  TO  2 
GO  TO  10 

19  PSII-PSII+TWOn 

TAVG-TCALC/TWOPI 

CALL  DWNWSH (V.TAVG.RHO, AREA, Al, ALPHA, W) 

TCALOO.  0 

IF(PSI.GE.PSILIM)GO  TO  18 
GO  TO  15 
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18  CONTINUE 

WRITE(*, *) ' PSILIM—' .PSILIM, '  KNTRL- ' , KNTRL 

IF(KNTRL.EQ.2)THEN 

WRITE(*,*) 'KNTRL-' .KNTRL 

PSII-PSIII 

PSI-PSII 

BETA-BETAI 

BETAD-BETADI 

W-WI 

TIME-TIMEI 
TC-TPRED 
FB-FB-SWTCH2 
GO  TO  21 
ENDIF 
KNTRL- 2 

IF(PSI.GE.PSIMAX)GO  TO  5000 
WRITE  (*,*)'  KNTRI^-2 ' 

20  CONTINUE 
PSII-PSI 
PSIPRI— PSIPRT 

PS I LIM-PS I +NPRED*TWOP I 

TIMEI-TIME 

BETAI-BETA 

BETADI-BETAD 

WI-W 

PSIII-PSI 
GO  TO  1 

8  IF(ABS(BETA) . GT . BETLIM) GO  TO  9 

SWTCH2-1.0 

GO  TO  5 

9  DELBT2-ABS (BETA) -BETLIM 

1 F ( DELBT1 . GE . DELBT2 ) GO  TO  14 

DELBT1-DELBT2 

GO  TO  5 

15  CONTINUE 

DELT1-0. 

TCALC-0 . 

GO  TO  6 

14  CONTINUE 

SWTCH2-0 . 0 
FB-FB+1 . 0 
SWITCH-1.0 

IF (COS (PSI) . GE . SIN (PS I) ) SWITCH-2 . 0 
WRITE(*,*) 'KNTRL-3' 

TC-TPRED* ( PS I - PS I I I ) /TWOPI/NPRED 

21  CONTINUE 
KNTRL-3 

PSI PRT-PS I PRI 

PSI-PSIII 

TIME-TIMEI 

BETA-BETAI 

BETAD-BETADI 

W-WI 

GO  TO  15 
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11  CONTINUE 

I F (TIME . LT . TPRINT ) GO  TO  12 

TH1DG-THETA1/DTR 

TH2DG— THETA2 /DTR 

REV— PSI/TWOPI 

BETADG-BETA/DTR 

FDBK1G-FDBK1/DTR 

FDBK2G-FDBK2 /DTR 

WRITE(*,*) 'TIME-' .TIME, '  REV-' ,REV,  '  BETA-' , BETADG 
WRITE (20,500 ) TIME , TH1DG , TH2DG , FDBK1G , FDBK2G , REV , BETADG 
IF(PSI.GE.PSIMAX)GO  TO  5000 
TPRINT-TPRINT+DELTPR 

12  CONTINUE 

IF (TIME . GE . TIMEI+TC ) GO  TO  13 
GO  TO  10 

13  CONTINUE 
KNTRL-2 
DELBT1-0 . 

GO  TO  20 

5000  CLOSE (UNIT-20) 

STOP 

END 

SUBROUTINE  SUBCNTRL 

COMMON/COM1/THE1LM , THE2LM , KNTRL, TCON , TOFF , THE1I , THE2I , SWITCH 
COMMON/ COM2/DELFB1 , DELFB2 , FB , PI , CRATE1 , CRATE2 , FDBK1 , FDBK2 
COMMON/COM5/PSI , EPS , CO , CT , THETAO , THETAT , THETAl , THETA2 , KBETA , BETA 
COMMON/COM6/FB1LIM, FB2LIM 
REAL  KBETA 

1  IF (TCON . GT . TOFF) GO  TO  2 
DELTH1-CRATE1*TC0N 
DELTH2-CRATE2*TCON 

IF(DELTH1 . GT . THE1LM) DELTH1-THE1LM 
IF(DELTH2 . GT . THE2LM)DELTH2-THE2LM 
THETA1-THE1I+DELTH1 
THETA2-THE2I +DELTH2 
GO  TO  3 

2  DELTH1-CRATE1* (TCON -TOFF) 

DELTH2-CRATE2* ( TCON - TOFF ) 

IF (DELTH1 . GT . THE1LM) DELTH1-THE1LM 
IF(DELTH2 . GT . THE2LM) DELTH2-THE2LM 
THETA1-THE1I+THE1LM- DELTH1 
THETA2-THE2I+THE2LM-DELTH2 

3  IF(FB.LE.O. )FB— 0.0 
FDBK1-FB*DELFB1 
FDBK2-FB*DELFB2 

IF(FDBK1 .GT. FB1LIM)FDBK1-F31LIM 
IF ( FDBK2 . GT . FB1LIM) FDBK2-FB2LIM 
IF( SWITCH . EQ . 1 . ) FDBK2-0 . 0 
I F( SWITCH . EQ . 2 . ) FDBK1-0 . 0 
THETA1-THETA1 - FDBK1 
THETA2-THETA2 - FDBK2 
IF(THETA1 . LT . THE1I )THETA1-THE1I 
IF(THETA2 . LT . THE2I )THETA2-THE2I 
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4  RETURN 

END 

SUBROUTINE  DWNWSH( V , TAVG , RHO , AREA , Al , ALPHA ,  W) 

W-SQRT(0 . 5*( -V**2+SQRT(V**4+(TAVG/RHO/AREA)**2) ) ) 

DO  1  1-1,10 

VPRIME-SQRT ( (V-W*SIN(ALPHA+Al) )**2+(W*COS(ALPHA+Al) )**2) 
W-TAVG/2/RHO/VPRIME/AREA 
1  CONTINUE 

RETURN 
END 

SUBROUTINE  SUBFLAP 

REAL  KBETA , MW . MIF , Ml , M2 , MAERO 

COMMON/COM3/CLI (20) , ALPHAL(20) , CDI(33) ,ALPHAD(33) ,ALPHAB,CL,CD 
C0MM0N/C0M4/0MEGA , V , ALPHA , W , BETAD , RHO , DELPSI , DELR , T , MIF , MW , D 
COMMON /COM5 /P  S I , EPS , CO , CT , THETAO , THETAT , THETA1 , THETA2 , KBETA , BETA 
PI-3.14159 
DELT-DELPS I/OMEGA 
S ALPHA-S IN ( ALPHA) 

CALPHA-COS (ALPHA) 

C  REM  START  OF  INTEGRATION  OF  THRUST  OVER  RADIUS  AT  A  PARTICULAR  PSI 

1  T-0 

SPSI— SIN(PSI) 

CPSI-COS(PSI) 

EPSR— EPS*D/2 
R-EPSR 
Ml— 0 
DTI— 0 
MAERO-O 

C  REM  RETURN  TO  HERE  AFTER  INCREMENTING  RADIUS,  R 

2  X— R/D*2 
C-CO- (C0-CT)*X 

C  THETA-BLADE  PITCH  ANGLE 

THETA-THETA0+THETAT*X+THETA1*CPSI+THETA2*SPSI+KBETA*BETA 
C  VTHETA-RESULTANT  TANGENTIAL  VELOCITY 

VTHETA-OMEGA*R+V*CALPHA*S  PSI 
C  VU-RESULTANT  VELOCITY  UP  THROUGH  DISK 

VU—V*S ALPHA - W- (R- EPSR) *BETAD - BETA*V*CALPHA*CPSI 
ALPHAB-ATAN ( ABS (VU/VTHETA) ) 

C  START  OF  LOGIC  TO  DETERMINE  CL  AND  CD  AND  RESOLUTION  FACTORS 
CLFAC-1 . 0 
CDFAC-1 . 0 

IF  ( VTHETA . GE . 0 . 0 . AND . VU . GE . 0 . 0 ) THEN 
ALPHAB-THETA+ALPHAB 
IF(ALPHAB . LT . 0 . 0)CLFAC— 1 . 

ALPHAB-ABS (ALPHAB) 

GO  TO  3 

ENDIF 

IF  (VTHETA. GT. 0.0. AND.VU.LT. 0.0)  THEN 
ALPHAB-THET A - ALPHAB 
IF (ALPHAB . LT . 0 . 0 ) CLFAC-- 1 . 

ALPHAB-ABS (ALPHAB) 

GO  TO  3 


ENDIF 
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IF  (VTHETA.LT. 0.0. AND.VU.LT. 0.0)  THEN 
ALPHAB-PI -ALPHAB -THETA 
CLFAC— 1.0 
CDFAC— 1.0 

IF (ALPHAB . GT . PI ) CLFAC-1 . 0 
IF (ALPHAB . GT . PI ) ALPHAB-2 . *PI -ALPHAB 
GO  TO  3 

ENDIF 

IF  ( VTHETA . LT . 0 . 0 . AND . VU . GT . 0 . 0 )  THEN 
ALPHAB— PI+THETA- ALPHAB 
CLFAC-1. 0 
CDFAC— 1.0 

IF(  ALPHAB .  GT .  PI)  CLFAC—  1 . 0 
I F( ALPHAB . GT . PI ) ALPHAB-2 . *PI -ALPHAB 

ENDIF 

C  END  OF  LOGIC  ON  CL,  CD  AND  RESOLUTION  FACTORS 

3  CALL  AIRFOIL 
CLr*CL*CLFAC 
CD-CD*CDFAC 
VRSQD-VU**2+VTHETA**2 
DLDR-RHO/2*VRSQD*CL*C 
DDDR— RHO/2*VRSQD*CD*C 

C  DLDR  AND  DDDR  ARE  LIFT  AND  DRAG  DERIVATIVES  WRT  TO  R 
VR-SQRT ( VRSQD) 

DT2-DLDR*ABS (VTHETA/VR) +DDDR*ABS (VU/VR) 

M2— (R-EPSR)*DT2 
DMDR— (M1+M2 ) /2 
DTDR— (DT1+DT2 ) /2 

C  DTDR  AND  DMDR  ARE  RADIAL  THRUST  AND  MOMENT  DERIVATIVES  WRT  R 
Ml— M2 
DT1-DT2 

MAERO-MAERO+DMDR*DELR 

T-T+DTDR*DELR 

C  T  AND  MAERO  ARE  BLADE  THRUST  AND  MOMENT  AT  INSTANT  OF  TIME 
R-R+DELR 

IF  (R.GT.D/2)  GO  TO  4 
GOTO  2 

4  BETADD-MAER0/MIF-BETA*0MEGA**2 -MW/MIF 
BETAD-BETAD+BETADD*DELT 
BETA-BETA+BETAD*DELT+BETADD*DELT**2/2 . 

Ml-0 . 

DT1-0 . 

PSI2-PSI 

PSI1-PSI2 

C  T  IS  THE  BLADE  THRUST  AT  AN  INSTANT  OF  TIME  AND  PSI 
RETURN 
END 

C  SUBROUTINE  FOR  CL  AND  CD  AS  FUNCTION  OF  ALPHA 
SUBROUTINE  AIRFOIL 

COMMON/COM3/CLI ( 20 ) ,ALPHAL(20) ,CDI(33) ,ALPHAD(33) , ALPHAB , CL , CD 
DO  1  1-2,20 

IF ( ALPHAB. EQ.ALPHAL(I) )  THEN 
CL-CLI(I) 


h  in  n  i*no  4 
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GO  TO  2 

ENDIF 

IF (ALPHAB . LE . ALPHAL( I ) )  THEN 
AI-CLI(I) 

BI-CLI(I-l) 

Cl— ALPHAL(I) 
DI-ALPHAL(I-l) 

GO  TO  5 

ENDIF 

CONTINUE 

CL-(AI-BI)*(ALPHAB-DI)/(CI-DI)+BI 

CONTINUE 

DO  3  1-2,33 

IF(ABS (ALPHAB) . EQ . ALPHAD ( I ) )  THEN 
CD-CDI(I) 

GO  TO  4 

ENDIF 

IF (ABS (ALPHAB) . LE . ALPHAD ( I ) )  THEN 
AI-CDI(I) 

BI— CDI(I-l) 

Cl— ALPHAD(I) 

DI— ALPHAD ( I - 1) 

GO  TO  6 

ENDIF 

CONTINUE 

CD- (AI - BI )* (ALPHAB - DI )/( Cl - DI ) +BI 

RETURN 

END 


FORTRAN  listing  for  PROGRAM  DCFLAP 
(digital  controller) 


Ci  o  o  o 


88 


PROGRAM  DCFLAP 

C - - - 

C  Program  to  simulate  the  helicopter  rotor  blade  flapping 

C  motion  with  real-time  control  system  preventing  excessive 

C  flapping.  The  control  system  can  be  cut  off  by  setting 

C  GAIN-0 . 

C  Flapping  angle  BETA  and  flapping  angular  velocity  BETAD 

C  as  a  function  of  NREV  (number  of  revolutions)  are  stored 

C  in  file  RESULT .  DAT . 

C  Pitch  angles  THETAO  (collective) ,  THETA1  (lateral  cyclic) , 

C  and  THETA2  (longitudinal  cyclic)  as  a  function  of  NREV  are 

C  stored  in  file  CONTROL.DAT. 

C  Feedback  gain  GAIN  of  the  control  system  as  a  function  of 

C  NREV  is  stored  in  GAIN.DAT.  (An  adaptive  gain  technique 

C  is  used.) 

C 

C  Aerospace  Engineering  Department 

C  Pennsylvania  State  University 

C  University  Park,  PA  16802 

C -  - -  - - - . - . . _____ - - - - 

c 

c 

CHARACTER  REPLY*1 , TRIM* 5 

REAL  KBETA,MW,MIF,M1 ,M2 ,MAERO,NMAX,NBET,MU, 

&  LAMDA , NREV , NPRED 

COMMON/PMT/PI , DTR , STEP , NMAX , PSIB , DPSIB , W 
COMMON/INPUT/THETO , THET1 , THET2 , ALPHA , BET1 , BET2 
COMMON/SIGNAL/TRIM . VAR1 , VAR2 , VAR3 , VAR4 
COMMON/PTB/PERTB , DO , D1 , D2 , DALPHA 

COMMON/S IMU/B , D , CO , CT , EPS , KBETA , THETAT , MH , MIF , V , VT , DELX 
& , RHO , KTEN , CRATEO , CRATE1 , CRATE2 , CRATEA , GAIN , SOS , BEEP ,U,U0, UPRED 
DATA  TRIM/' EXIST'/ 

C 

C -  —  . .  ■■  - - - - - - ________ 

C  FILE  ' RESULT . DAT '  STORES  THE  VALUES  OF  PS IB, BETA, BETA  DOT  (IN  DEGREES) 
C  FILE  ' TRIM . DAT '  STORES  THE  CONTROL  INPUTS  FOR  THE  CURRENT  TRIM 
C  FILE  ' IC.DAT'  STORES  THE  INITIAL  CONDITION  TABLE 
C 

OPEN (100 , FILE- ’ RESULT . DAT ' , STATUS- ’ NEW' ) 

OPEN( 101 , FILE- ' TRIM . DAT ' , STATUS- ' UNKNOWN ' ) 

OPEN(102 , FILE- ’IC.DAT', STATUS- ' UNKNOWN' ) 

OPEN ( 107 , FILE- ' GAIN . DAT ' , STATUS- ' UNKNOWN ' ) 

REWIND  3 
REWIND  101 
REWIND  102 


NBET  ---  THE  NUMBER  OF  POSITIONS  PER  REVOLUTION  TO  BE  CALCULATED 
NBET-72 

STEP  --  STEP  IN  THE  INITIAL  CONDITION  TABLE 
STEP-360. /NBET 

C  - 

PI  -  3.14159265 
DTR-PI/180. 
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C  - — .  . - . . . . . — . . — — - - 

C  HAVE  YOU  GENERATED  THE  TRIM  (OR  NORMAL)  FLIGHT  CONDITIONS? 

C  NO- ->GOTO  1010  TO  INPUT  THE  INITIAL  VALUES  NEEDED  TO  GENERATE: 

C  ' TRIM . DAT '  &  ' IC . DAT ' 

C  YES- ->READ  THE  INITIAL  VALUES  FROM  THE  TWO  FILES  ' TRIM . DAT '  &  'IC.DAT' 
C  THEN,  GOTO  1020 

C 

1000  WRITE(6,*)  'Have  you  generated  the  trim  condition? (Y/N) : ' 

READ(5 , 31)  REPLY 

IF(REPLY.EQ. 'N' . OR . REPLY . EQ . 'n')  THEN 
TRIM- 'NONE' 

GOTO  1010 

ELSE  IF(REPLY.EQ. 'Y' . OR . REPLY . EQ . ' y')  THEN 
READ (101,32)  THETOD , THET1D , THET2D , ALPHAD , W 
32  FORMAT ( IX , 5G15 . 8 ) 

READ(102 , 35)  PSIBD , BET10D , BET20D 
REWIND  101 
REWIND  102 
GOTO  1020 
ELSE 

GOTO  1000 
ENDIF 

1010  WRITE(6,*)  'Enter  one  blade  azimuth  angle  in  degrees:' 

READ(5 ,*)  PSIBD 

WRITE(6,*)  'Enter  collective  in  degrees:' 

READ(5 ,*)  THETOD 

WRITE(6 ,*)  'Enter  longitudinal  cyclic  in  degrees:' 

READ(5 ,*)  THET2D 

WRITE(6 ,*)  'Enter  lateral  cyclic  in  degrees:' 

READ (5,*)  THET1D 

WRITE(6 ,*)  'Enter  ALPHA  in  degrees:' 

READ (5 , *)  ALPHAD 
BET10D— 0 . 

BET20D-0 . 

W-0.0 

GOTO  1300 

C  -  -  - - -  ,  —  M—  .  - - - - - - - - - - 

C  THE  FOLLOWING  BLOCK  ALLOWS  YOU  TO  EXAMINE  THE  EFFECTS  OF  EACH  CONTROL 
C  INPUTS  OR  DISTURBANCES  ON  THE  ROTOR  FLAPPING 
C 

1020  WRITE(6,*)  'These  are  the  current  inputs:' 

WRITE(6 , 34)  THETOD, THET2D , TH.ET1D , ALPHAD 
34  FORMAT( IX, 'Collective  is ' ,G15 . 8 ,' degrees ' 

&  /IX, 'Longitudinal  cyclic  is' ,G15. 8, 'degrees' 

&  /IX, 'Lateral  cyclic  is' ,G15. 8, 'degrees' 

&  /IX, 'ALPHA  is' ,G1S. 8, 'degrees') 

IF(TRIM.EQ. 'NONE')  THEN 
PERTB-1 . E+30 
GOTO  1205 
ELSE 

1100  WRITE(6,*)  'Do  you  want  to  give  perturbations? (Y/N) : ' 

READ(5 , 31)  REPLY 

IF(REPLY.EQ. 'Y' . OR . REPLY . EQ . 'y')  THEN 
GOTO  1101 
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ELSE  IF(REPLY.EQ. 'N' .OR. REPLY. EQ. 'n')  THEN 

GOTO  1300 

ELSE 

GOTO  1100 
END  IF 

1101  WRITE(6,*)  'Enter  number  of  revolutions  at  which  the  ', 

&' perturbations  are  given:' 

READ(5 ,*)  PERTB 

1200  WRITE(6,*)  'Enter  perturbations  on  Collective,  Longitudinal', 

&' Cyclic,  Lateral  Cyclic,  and  ALPHA  in  degrees:' 

READ(5 ,*)  DO , D1 , D2 , DALPHA 
D0-D0*DTR 
Dl-D1*DTR 
D2-D2*DTR 
DALPHA-DALPHA*DTR 
31  FORMAT (Al) 

C . 

WRITE(6,*)  'Enter  changing  rates  on  Collective,  Longitudinal', 

&' Cyclic,  Lateral  Cyclic,  and  ALPHA  in  deg/sec:' 

READ ( 5 , * )  CRATEO , CRATE1 , CRATE2 , CRATEA 
C  Transform  CRATEO , CRATE1 , CRATE2 , CRATEA  to  RAD/SEC 
CRATE0-CRATE0*DTR 
CRATE1-CRATE1*DTR 
CRATE2-CRATE2*DTR 
CRATEA-CRATEA*DTR 
C 

WRITE (6,*)  'INPUT  FEEDBACK  GAIN:' 

READ(5 ,*)  GAIN 
ENDIF 

1205  CONTINUE 

C . 

c 

c  —  ■  —  ■■  . .  -  ..  . . . . — 

C  THE  FOLLOWING  BLOCK  PERFORMS  LINEAR  INTERPOTATION  IN  THE  INITIAL  CONDITION 
C  TABLE  TO  GET  THE  CORRECT  INITIAL  CONDITION  FOR  THE  PSIB  INPUT 
C 

C  NENTRY--  THE  ENTRY  PLACE  NUMBER- 1  IN  THE  INITIAL  CONDITION  TABLE 

1206  NENTRY«*JNINT(PSIBD/STEP-  .5) 

DO  1310  JUMP-1, NENTRY 

1310  READ(102 , 35)  VNUL 

READ(102 , 35)  PSIBD1,BET1D1,BET2D1 
READ(102 , 35)  PSIBD2 , BET1D2 , BET2D2 

BET10D-BET1D1+ ( PS IBD - PS IBD1 ) / ( PS I BD2 - PSIBD1 ) * ( BET1D2 - BET1D1 ) 
BET20D-BET2D1+(PSIBD-PSIBD1)/(PSIBD2-PSIBD1)*(BET2D2-BET2D1) 

35  F0RMAT(1X, 3G15 . 8) 

REWIND  102 
C 
C 

C  ■  ■  — - - - - - - 

1300  PSIB-PSIBD*DTR 

THETO-THETOD*DTR 
THET1-THET1D*DTR 
THET2 -THET2 D*DTR 
BET10-BET10D*DTR 
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BET20-BET20D*DTR 

ALPHA-ALPHAD*DTR 

NMAX-10. 

DPS IB-2. *PI/NBET 
IF(TRIM.EQ. 'EXIST')  THEN 
NREV-PSIB/2./PI 

WRITE(100 , 35)  NREV , BET10D , BET20D 
ENDIF 
C 

CALL  BWFLAP 
REWIND  102 
C 

C  THE  FOLLOWING  WRITE  STATEMENT  RECORDS  THE  CONTROL  INPUTS  FOR  THE  CURRENT 
C  TRIM  CONDITION 

IF (TRIM. EQ. 'NONE' )  THEN 

WRITE( 101 , 32 )  THETOD , THET1D , THET2D , ALPHAD , W 
REWIND  101 
WRITE (6 , 36) 

36  FORMAT ( 30 (/) ,10X, 'TRIM  GENERATED! ', 5 (/) ) 

TRIM-' EXIST' 

GOTO  1020 
ENDIF 
STOP 
END 

C - - - - - - - - - - - - -  - . . — 

SUBROUTINE  WRITE 
REAL  NMAX, NREV 
CHARACTER  TRIM*5 

COMMON /PMT /PI , DTR, STEP , NMAX, PSI , DELPSI , W 
COMMON/INPUT/THETAO , THETAl , THETA2 , ALPHA , BETA , BETAD 
COMMON/SIGNAL/TRIM , NREV , VAR1 , VAR2 , VAR3 
PSID— PSI/DTR 
BET-BETA/DTR 
BETD-BETAD/DTR 
IF(TRIM.EQ. 'EXIST')  THEN 
WRITE (100, 35)  NREV , E  ET , BETD 
35  FORMAT ( IX , 3G15 . 8 ) 

ELSE  IF(PSI.GE. (16.*PI-STEP*PI/360.) .AND.PSI.LE. (18.*PI))  THEN 

PSI1ST-(PSI-16.*PI)*180./PI 

WRITE (102 ,35)  PSI1ST, BET, BETD 

ENDIF 

RETURN 

END 


C 

SUBROUTINE  BWFLAP 

CHARACTER  BEEP*1 , DETECT* 1 , TRIM*5 

REAL  UTHET1 ( 2 ) , UTHET2 ( 2 ) 

REAL  KBETA , MW , MIF , Ml , M2 , MAERO , NMAX , NBET , MU , 
&  LAMDA , NREV , NPRED 
DIMENSION  U(10),U0(10) 

COMMON/PMT/PI , DTR , STEP , NMAX , PS I , DELPS I , W 
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COMMON/INPUT/THETAO , THETA1 , THETA2 , ALPHA , BETA , BETAD 
COMMON  CLI(20) ,ALPHL(20) ,CDI(33) ,ALPHD(33) , ALPHAB , CL , CD 
COMMON/SIMU/B , D , CO , CT , EPS , KB ETA , THETAT , MW , MIF , V , VT , DELX 
& , RHO , KTEN , CRATEO , CRATE 1 , CRATE2 , CRATEA , GAIN , SOS , BEEP , U , UO , NPRED 
COMMON/PTB/PERTB , DO , D1 , D2 , DALPHA 
COMMON /HOLD/UTHET 1 .UTHET2 
COMMON/SIGNAL/TRIM , NREV , A1 , B1 , AREA 

COMMON/CLOSED/T1 , T2 , T3 , T4 , FI , F2 , F3 , F4 , DENOM , All , Al2 , Al 3 , A14 , 

&  Bll.TAU.CBAR.GAMMAF, SIGMA, ASIG 
C  NPRED  is  number  of  predict-ahead  revolutions 
NPRED- 3 . 


C 


1 


3 


5 


6 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DETECT-' 0' 

OPEN (UNIT-1 , FILE- ' CL . DAT ' , STATUS- ' OLD ' ) 
DO  1  1-1,20 

READ(1,*)  ALPHL(I) ,CLI(I) 

CONTINUE 
CLOSE  (UNIT-1) 

OPEN (UNIT-2 , FILE- ' CD . DAT ' , STATUS- ' OLD ' ) 
DO  3  1-1,33 

READ(2 ,*)  ALPHD(I) ,CDI(I) 

CONTINUE 
CLOSE (UNIT-2) 

PI-3.14159 

DO  5  1-1,20 

ALPHL ( I ) -ALPHL ( I ) *DTR 

CONTINUE 

DO  6  1-1,33 

ALPHD ( I ) -ALPHD ( I ) *DTR 

CONTINUE 

PSII-0.0 

PROGRAM  READS  IN  ORDER: 


FLP5.DAT 

1.  NUMBER  OF  BLADES 

2.  ROTOR  DIAMETER 

3.  ROOT  CHORD  AT  R-0 

4.  TIP  CHORD 

5.  DIMENSIONLESS  HINGE  OFFSET 


FLP6 . DAT 

1.  PITCH-FLAP  COUPLING,  KBETA 

2.  TOTAL  TWIST,  DEGS . .NEGATIVE  FOR  WASHOUT 

3.  BLADE  WEIGHT  MOMENT  ABOUT  FLAPPING  AXIS 

4.  MASS  MOMENT  OF  INERTIA  ABOUT  FLAPPING  AXIS 

5.  TIP  SPEED  IN  FPS 


FLP7 . DAT 

1.  DIMENSIONLESS  INCREMENT  IN  RADIUS,  DELX 

2.  INCREMENT  IN  AZIMUTH  ANGLE, DEGS,  DELPSI 

3.  NUMBER  OF  DELPSI  BETWEEN  PRINTOUTS 


FLP8.DAT 

1.  FORWARD  SPEED  IN  KTS . 

2.  AIR  DENSITY,  SLUGS/CU.FT. 

3.  COLLECTIVE  PITCH,  THETAO , DEGREES 


93 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


4.  LATERAL  PITCH , THETAl ,  DEGREES 

5.  LONGITUDINAL  PITCH , THETA2 . DEGREES 

6.  DISC  PLANE  ANGLE  OF  ATTACK, ALPHA, DEG S . 

FLP5.DAT  2,44,2.25,2.25, .01 

FLP6.DAT  0,-10,3122,1422,738 

FLP7.DAT  .05,5,5 

FLP8.DAT  80, .002378,15.1,1.9, -1.66, -4.63  (FWD  FLIGHT) 

FLP8.DAT  0, .002378,15,0,0,0  (HOVER) 


OPEN (UNIT-5 , FILE- 'FLP5.DAT', STATUS- ' OLD ' ) 

READ  (5,*)B,D,C0,CT,EPS 
CLOSE (UNIT-5) 

OPEN (UNIT-6 , FILE- ' FLP6 . DAT ' , STATUS- ' OLD ' ) 

READ  (6 ,*)KBETA,THETAT,MW,MIF,VT 
CLOSE(UNIT-6) 

OPEN (UNIT-7 , FILE- ' FLP7 . DAT ' , STATUS- ' OLD ' ) 
READ  (7,*)  DELX 
CLOSE (UNIT-7) 

OPEN (UNIT-8 , FILE- ' FLP8.DAT' , STATUS- ' OLD ' ) 
READ  (8,*)  V.RHO 
CLOSE (UNIT-8) 

OPEN ( 105 , FILE- ' CONTROL . DAT ' , STATUS- ' UNKNOWN ’ ) 
V-V*l . 69 
AREA— P I *D*D/4 . 

OMEG A— VT /D*  2 . 

DELT-DELPS I /OMEGA 
DELR-DELX*D/2 . 

PSIO-PSI 
TAVG— 0 . 

TTAVG-0 . 

C  KTEN  specifies  the  sampling  period 
KTEN-30 . +DTR/DELPSI+0 . 5 
C 

IKTEN-0 

THETAT-THETAT*DTR 

MU-V/VT 

Tl-(.941+MU**2/2)/2 
T2-( . 941/3 . +MU**2/2 . )* . 97 
T3- . 941/4 . * ( . 941+MU**2 ) 

T4-MU/2 . * ( . 941+MU**2/4 . ) 

F1-. 97**3/3 . 

F2- . 941/4 . * ( . 941+MU**2 ) 

F3- . 97**3*( . 941/5 . +MU**2/6 . ) 

F4-MU* . 97**3/3 . 

DENOM- . 941 -MU**2/2 . 

All-4 .*(MU* . 941/2 . -MU**3/8 . )/ . 941/DENOM 

A12-8 . *MU* . 97/3 . /DENOM 

A13-2 . *MU* . 941/DENOM 

A14— ( . 941+3 ./2 .*MU**2) /DENOM 

Bll-4 . *MU* . 97/3 . /( . 941+MU**2/2 . ) 

TAU-MW /MI F/OMEG A** 2 
CBAR- ( CO+CT ) /2 . 

GAMMAF-CBAR*RHO*  5 . 7*D**4/32 . /MIF 


OOU  VOOCSIUUOOOUOO  u  oo 
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SIGMA-B*CBAR/PI/D*2 . 
ASIG— 5 . 7*SIGMA 
BETAO-O . 0 
Bl-0.0 


WRITE  INITIAL  CYCLIC  INPUTS  AND  ZERO  CONTROL  INPUTS  TO  FILE  "CONTROL. DAT" 

NREV-PSI/2./PI 
CNTRL1-0 . 

CNTRL2-0 . 

E1-THETA1/DTR 

E2-THETA2/DTR 

WRITE ( 105 , 6 1 )  NREV , CNTRL1 , CNTRL2 , El , E2 
1  FORMAT ( IX , 5G15 . 8 ) 


THRST1-0 . 
280  CONTINUE 


START  OF  INTEGRATION  LOOP  FOR  PS I 
RETURN  TO  HERE  AFTER  INCREMENTING  PSI 


THIS  BLOCK  STIPULATES  THE  PILOT  INPUTS. 

CRATEO , CRATE1 , CRATE2 , CRATEA  ARE  THE  INPUT  RATES  OF  THETAO , THETAl , THETA2 , ALPHA 
RESPECTIVELY 

NREV— PS 1/2 . /PI 
IF (NREV . GE . PERTB)  THEN 
IF(DETECT.EQ. '0' )  THEN 
THETAO-THETAO+DO 
THETAl-THETAl+D 1 
THETA2-THETA2+D2 
ALPHA-ALPHA+DALPHA 
DETECT-' 1’ 

ELSE 

THETAO-THETAO+CRATEO*DELPS I/OMEGA 
THETAl-THETAl+CRATE 1*DELP S I /OMEGA 
THETA2-THETA2+CRATE2*DELPS I/OMEGA 
ALPHA-ALPHA+CRATEA*DELP  S I /OMEGA 
ENDIF 
ENDIF 

UTHET1 ( 1 ) -UTHET1 ( 2 ) 

UTHET1 ( 2 ) -THETAl 
UTHET2 ( 1 ) -UTHET2 ( 2 ) 

UTHET2 ( 2 ) -THETA2 


SPSI-SIN(PSI) 

CPSI-COS(PSI) 

SALPHA-SIN (ALPHA) 

CALPHA-COS (ALPHA) 

LAMDA-MU*ALPHA 
DO  9050  1-1,20 

CTR-ASIG/2 . * ( LAMDA*T1+ (THETA0+KBETA*BETA0 ) *T2+THETAT*T3+ ( THETA2 
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1-KBETA*B1)*T4) 

C  FAC— V**4+ ( TTAVG /RHO /AREA ) **  2+W**  3*2.*V*SIN( Al+ALPHA ) 

IF(CFAC . LT . 0 . )  THEN 

W-SQRT(0 . 5*( -V**2+SQRT(V**4+(TTAVG/RHO/AREA)**2) ) ) 

GO  TO  9051 
END  IF 

WFAC— V**2+SQRT ( CFAC ) 

IF(WFAC.LT.O.)THEN 

W-SQRT(0 . 5*( -V**2+SQRT(V**4+(TTAVG/RHO/AREA)**2) ) ) 

GO  TO  9051 
END  IF 

W-SQRT (0 . 5* ( - V**2+SQRT (CFAC) ) ) 

9051  WVT-W/VT 

LAMDA-MU* ALPHA- WVT 

BETAO-GAMMAF* ( LAMDA*F1+ ( THETA0+KBETA*BETA0 ) *F2+THETAT*F3+ ( THETA2 - 
1KBETA*B1)*F4) -TAU 

Al-LAMDA*Al 1+ ( THETAO+KBETA*B  ETAO ) * Al2+THETAT*Al 3+ ( THETA2 - KBETA*B 1 
1 ) *Al4 

B1-BETA0*B1 1 - ( THETA1 - KBETA*A1 ) 

9050  TTAVG-RHO*AREA*VT**2*CTR 
THRST2-0 

C  START  OF  INTEGRATION  OVER  X  FROM  XH  TO  1 
EPSR— EPS*D/2 
R-EPSR 
Ml-0 
DT1-0 
MAERO-O 

C  RETURN  TO  HERE  AFTER  INCREMENTING  RADIUS,  R 
2380  X— R/D*2 

C— CO- (C0-CT)*X 

THETA-THETA0+THETAT*X+THETAl*CPSI+THETA2*SPSI+KBETA*BETA 
VTHET A-OMEGA*R+ V*CALPHA*  S  P  S I 

VU— V*SALPHA - W - (R- EPSR) *BETAD - BETA*V*CALPHA*CPSI 
ALPHAB-ATAN (ABS (VU/VTHETA) ) 

CLFAC-1 . 0 
CDFAC-1 . 0 

IF  ( VTHETA . GE . 0 . 0 . AND . VU . GE . 0 . 0 ) THEN 
ALPHAB-THETA+ALPHAB 
IF(ALPHAB .  LT .  0 . 0 )  CLFAC—  1 . 

ALPHAB-ABS ( ALPHAB ) 

GO  TO  2381 

ENDIF 

IF  (VTHETA. GT. 0.0. AND. VU.LT. 0.0)  THt., 

ALPHAB-THETA - ALPHAB 
IF  (ALPHAB .  LT .  0 . 0  )  CLFAC— 1 . 

ALPHAB-ABS (ALPHAB) 

GO  TO  2381 

ENDIF 

IF  (VTHETA. LT. 0.0. AND. VU.LT. 0.0)  THEN 
ALPHAB-PI -ALPHAB-THETA 
CLFAC— 1.0 
CDFAC— 1.0 

IF (ALPHAB . GT . PI ) CLFAC-1 . 0 
IF(ALPHAB . GT . PI ) ALPHAB-2 . *PI -ALPHAB 
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2381 


2690 


C 

C 


GO  TO  2381 

ENDIF 

IF  ( VTHETA . LT . 0 . 0 . AND . VU . GT . 0 . 0 )  THEN 
ALPHAB-PI+THETA-ALPHAB 
CLFAC-+1 . 0 
CDFAC--1.0 

IF(ALPHAB .  GT .  PI )  CLFAC—  1 . 0 

IF ( ALPHAB . GT . PI ) ALPHAB-2 . *PI - ALPHAB 

ENDIF 

CALL  AIRFOIL 

CL-CL*CLFAC 

CD— CD*CDFAC 

VRSQD— VU**2+VTHETA**2 

DLDR-RHO/2*VRSQD*CL*C 

DDDR-RHO/2*VRSQD*CD*C 

VR-SQRT ( VRSQD) 

DT2-DLDR*ABS ( VTHETA/VR ) +DDDR*ABS (VU/VR) 

M2— (R-EPSR)*DT2 
DMDR— (M1+M2 ) /2 
DTDR— (DT1+DT2 ) /2 
Ml— M2 
DT1-DT2 

MAERO-MAERO+DMDR*DELR 

THRST2-THRST2+DTDR*DELR 

R-R+DELR 

IF  (R.GT.D/2)  GO  TO  2690 
GOTO  2380 

BETDD2-MAERO/MIF-BETA*OMEGA**2 -MW/MIF 
BETAD2-BETAD2+ ( BETDD1+BETDD2 )/2. *DELT 
BETAD-BETAD2 

BETA-BETA+ ( BETAD1+BETAD2 ) /2 . *DELT 

BETDD1-BETDD2 

BETAD1-BETAD2 

Ml-0 . 

DTI— 0 . 

PSI2-PSI 

PSI1-PSI2 

THRUST-  (  THRST 1+THRS  T  2 )  /  2  . 

THRST1-THRST2 
TAVG-TAVG+THRUST*DELPS I 

CALL  WRITE 

PSI-PSI+DELPSI 
PSII-PSII+  DELPSI 
IF  (PSII.LE.2.*PI)  GO  TO  9801 
TAVG— TAVG/2 . /PI*B 
CTRW-TAVG/RHO/AREA/VT** 2 
DO  9800  1-1,5 

CFAC-V**4+ ( TAVG/RHO/AREA) **2+W**3*2 . *V*SIN(Al+ALPHA) 
IF(CFAC.LT.O. )THEN 

W-SQRT ( 0 . 5* ( - V**2+SQRT (V**4+ (TAVG/RHO/AREA) **2 ) ) ) 
GO  TO  9800 
ENDIF 
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WFAC- - V**  2+SQRT ( CFAC ) 

IF ( WFAC .  LT .  0 . ) THEN 

W-SQRT ( 0 . 5*( -V**2+SQRT(V**4+(TAVG/RH0/AREA)**2) ) ) 

GO  TO  9800 
ENDIF 

W-SQRT ( 0 . 5* ( - V**2+SQRT ( CFAC ) ) ) 

9800  WVT-W/VT 
TAVG-0 . 

PSII-0. 

9801  IF (PS I . GT.NMAX*2 .*PI)  GOTO  100 
IKTEN-IKTEN+1 

IF(IKTEN. EQ.KTEN)  THEN 
C 

IF(TRIM.EQ. 'EXIST' .AND. ABS(GAIN) .GT.l.E-10)  CALL  CNTRL 
C 

IKTEN-0 
ENDIF 
GOTO  2280 
100  RETURN 
END 
C 
C 

C  SUBROUTINE  FOR  CL  AND  CD  AS  FUNCTION  OF  ALPHA 
SUBROUTINE  AIRFOIL 

COMMON  CLI ( 20 ) , ALPHL (20),CDI(33), ALPHD (33), ALPHAB , CL , CD 
DO  3040  1-2,20 
IF (ALPHAB . EQ . ALPHL ( I ) )  THEN 
Cb-CLI(I) 

GO  TO  4000 

ENDIF 

IF(ABS (ALPHAB) .LE.ALPHL(I))  THEN 
AI-CLI(I) 

BI-CLI(I-l) 

CI-ALPHL(I) 

DI-ALPHL(I-l) 

GO  TO  3050 

ENDIF 

3040  CONTINUE 

3050  CL-(AI -BI)* (ALPHAB -DI)/ (Cl -DI)+BI 
CL-CL*ABS (ALPHAB ) /ALPHAB 
4000  CONTINUE 

DO  3100  1-2,33 

IF (ABS (ALPHAB ) . EQ. ALPHD(I) )  THEN 
CD-CDI(I) 

GO  TO  3200 

ENDIF 

IF( ABS (ALPHAB) . LE.ALPHD(I) )  THEN 
AI-CDI(I) 

BI-CDI(I-l) 

CI-ALPHD(I) 

DI-ALPHD(I-l) 

GO  TO  3110 

ENDIF 

3100  CONTINUE 
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3110  CD— (AI-BI)*(ALPHAB-DI)/(CI-DI)+BI 
3200  RETURN 
END 
C 
C 

SUBROUTINE  CNTRL 

CHARACTER  BEEP*1 ,TRIM*5 

REAL  U(10) ,U0(10) ,UTHET1(2) ,UTHET2(2) 

REAL  KBETA,MW,MIF,M1 ,M2 , MAERO , NMAX , NBET , 

&  LAMDA , NREV , NPRED 
COMMON/PMT/PI , DTR , STEP , NMAX , PS 10 , DELPSI , WO 
COMMON/INPUT /THETAO , THETAl , THETA2 , ALPHA , XBETAO , XBETADO 
COMMON/SIMU/B , D , CO , CT , EPS , KBETA , THETAT , MW , MIF , V , VT , DELX 
& , RHO , KTEN , CRATEO , CRATE1 , CRATE2 , CRATEA , GAIN , SOS , BEEP , U , UO , NPRED 
COMMON  CLI(20) ,ALPHL(20) ,CDI(33) ,ALPHD(33) , ALPHAB , CL , CD 
COMMON/SIGNAL/TRIM , NREV , XAl , XB 1 , AREA 
COMMON/HOLD /UTHET1 .UTHET2 

C0MM0N/CL0SED/T1 , T2 , T3 , T4 , FI , F2 , F3 , F4 , DENOM , All , A12 , A13 , A14 , 

&  B11,TAU,CBAR,GAMMAF, SIGMA, ASIG 
Al— XAl 
B1-XB1 

BETAO-XBETAO 

BETADO-XBETADO 

E1-THETA1/DTR 

E2-THETA2/DTR 

W-WO 

BETA-BETAO 
BETAD-BETADO 
S ALPHA— S IN ( ALPHA ) 

CALPHA-COS (ALPHA) 

OMEGA— VT/D*2 

DELT-DELPSI/OMEGA 

TAVG-0 

DELR-DELX*D/2 

Ct  t  ++++++++++++++++++++++++++4  MM  H  +++++++  H  <  f -H-H-H-H-H 
WNREI^SQRT ( 1 . +1 . 5*EPS/ ( 1 . - EPS ) ) 

DAMP-. 5*(C0+CT)*RH0*5 . 73* (D/2 . )**4/16 ,/MIF*(l . -EPS)**4* 
&(1.+1./3.*EPS)/(1.-EPS)/WNREL 

CDELAY- (WNREL**2 - 1 . ) /SQRT ( ( WNREL**2 - 1 . ) **2+4 . * ( DAMP*WNREL) **2 ) 
DELAY-PI/2 . -ACOS (CDELAY) 

C  i  i  h-h  i  i  i  i  1 1 1 1  h  M  H  i  n  -H  i  m  1 1 1  1 1  i  i  1 1  i  i  h  i  n  i  h  1 1  m  1 1  n 
PS I00-PS 10 -KTEN*DELPS I 
PSI-PSI00 
OVERO-O . 0 
THRST1-0 . 

2280  CONTINUE 

SPSI-SIN(PSI) 

CPSI-COS(PSI) 

SALPHA-SIN(ALPHA) 

CALPHA-COS (ALPHA) 

C . 

C  This  block  is  the  First-Order  Hold 

THETAl-UTHET 1 ( 2 ) + (UTHET1 ( 2 ) -UTHET1 ( 1 ) ) /DELPSI* 

&  (KTEN*DELPS 1/ (NPRED*2 .*PI))*(PSI-PSI00) 
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THETA2-UTHET2 ( 2 ) + (UTHET2 ( 2 ) -UTHET2 ( 1 ) ) /DELPS I* 

&  (KTEN*DELPSI/(NPRED*2.*PI))*(PSI-PSI00) 

C . 

c 

LAMDA-MU*ALPHA 
DO  9050  1-1,20 

CTR-AS IG/2 .*( LAMDA*T 1+ ( THETA0+KBETA*BETA0 ) *T2+THETAT*T3+ ( THETA2 
1-KBETA*B1)*T4) 

CFAC— V**4+(TTAVG/RHO/AREA)**2+W**3*2 . *V*SIN (Al+ALPHA) 
IF(CFAC.LT.O. )  THEN 

W-SQRT ( 0 . 5*( - V**2+SQRT (V**4+(TTAVG/RH0/AREA) **2 ) ) ) 

GO  TO  9051 
END  IF 

WFAC— - V**2+SQRT ( CFAC ) 

IF(WFAC.LT.O. )  THEN 

W-SQRT (0 . 5*( -V**2+SQRT(V**4+(TTAVG/RHO/AREA)**2) ) ) 

GO  TO  9051 
END  IF 

W-SQRT ( 0 . 5* ( - V**2+SQRT ( CFAC ) ) ) 

9051  WVT-W/VT 

LAMDA— MU*ALPHA - WVT 

BETAO-GAMMAF* (  LAMDA*F1+ ( THETAO+KBETA*BETAO ) *F2+THETAT*F3+ ( THETA2 - 
1KBETA*B1 ) *F4 ) - TAU 

Al-LAMDA*A11+ ( THETAO+KBETA*BETAO ) *A1 2+THETAT*Al 3+ ( THETA2 - KBETA*B1 
1 ) *A14 

Bl— BETA0*B11 - (THETAl -KBETA*A1 ) 

9050  TTAVG-RHO*AREA*VT**2*CTR 
THRST2-0 

C  START  OF  INTEGRATION  OVER  X  FROM  XH  TO  1 
EPSR— EPS*D/2 
R-EPSR 
Ml— 0 
DTI— 0 
MAERO-O 

C  RETURN  TO  HERE  AFTER  INCREMENTING  RADIUS,  R 
2380  X-R/D*2 

C-CO- (C0-CT)*X 

THETA-THETA0+THETAT*X+THETA1*CPSI+THETA2*SPSI+KBETA*BETA 
VTHETA— OMEG A*R+V *  C ALPHA*  S  P  S I 

VU-V*SALPHA-W- (R- EPSR) *BETAD - BETA*V*CALPHA*CPSI 
ALPHAB-ATAN ( ABS (VU/VTHETA) ) 

CLFAC-1 . 0 
CDFAC-1 . 0 

IF  (VTHETA. GE. 0.0. AND. VU.GE.0.0)THEN 
ALPHAB-THETA+ALPHAB 
I F  ( ALPHAB .  LT .  0 . 0  )  CLFAC— 1 . 

ALPHAB-ABS (ALPHAB) 

GO  TO  2381 

END  IF 

IF  (VTHETA. GT. 0.0. AND.VU.LT. 0.0)  THEN 
ALPHAB-THETA - ALPHAB 
IF (ALPHAB . LT . 0 . 0)CLFAC— 1 . 

ALPHAB-ABS (ALPHAB) 

GO  TO  2381 


ENDIF 

IF  (VTHETA.LT. 0.0. AND.VU.LT. 0.0)  THEN 
ALPHAB-P I - ALPHAB - THETA 
CLFAC— 1.0 
CDFAC— 1.0 

IF (ALPHAB . GT . PI) CLFAC-1 . 0 
IF (ALPHAB . GT . PI ) ALPHAB-2 . *PI - ALPHAB 
GO  TO  2381 

ENDIF 

IF  ( VTHETA .  LT .  0 . 0 .  MD .  VU .  GT .  0 . 0 )  THEN 
ALPHAB-P I +TH ETA - ALPHAB 
CLFAC-+1 . 0 
CDFAC— 1.0 

IF  (ALPHAB .  GT .  PI )  CLFAC— 1 . 0 
IF (ALPHAB . GT . PI ) ALPHAB-2 . *PI -ALPHAB 

ENDIF 

2381  CALL  AIRFOIL 
CL-CL*CLFAC 
CD— CD*CDFAC 
VRSQD— VU**2+VTHETA**2 
DLDR— RHO/2*VRSQD*CL*C 
DDDR— RHO/2*VRSQD*CD*C 
VR-SQRT (VRSQD) 

DT2-DLDR*ABS (VTHETA/VR) +DDDR*ABS ( VU/VR) 

M2— (R-EPSR)*DT2 
DMDR— (M1+M2 ) /2 
DTDR— (DT1+DT2 ) /2 
Ml— M2 
DT1-DT2 

MAERO-MAERO+DMDR*DELR 

THRST2-THRST2+DTDR*DELR 

R-R+DELR 

IF  (R.GT.D/2)  GO  TO  2690 
GOTO  2380 

2690  BETDD2-MAERO/MIF-BETA*OMEGA**2 -MW/MIF 
BETAD2-BETAD2+(BETDD1+BETDD2)/2.*DELT 
BETAD-BETAD2 

BETA-BETA+ ( BETAD1+BETAD2 )/2. *DELT 

BETDD1-BETDD2 

BETAD1-BETAD2 

Ml-0. 

DT1-0 . 

PSI2-PSI 

PSI1-PSI2 

THRUST-(THRSTl+THRST2)/2 . 

THRST1-THRST2 
TAVG-TAVG+THRUST*DELPS I 
C 

BETLIM— l.*DTR 
IF (BETA . LT . BETLIM)  THEN 
OVER-BETA- BETLIM 
IF(ABS(OVER) .LT.ABS(OVERO) )  THEN 
C  SOS  IS  THE  ANGLE  AT  WHICH  MAXIMUM  OVERFLAPPING  OCCURS 

SOS-PS I+DELAY - DELPS I - AINT ( ( PSI+DELAY-DELPSI ) /2 . /PI ) *2 . *PI 
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BEEP-'Y' 

U(10)— ABS (OVERO) 

DO  210  IU-1,9 
210  U(IU)-U0(IU+1) 

CALL  FDBACK( OVERO) 

CNTRL1- - GAIN*OVERO*S IN ( SOS ) /DTR 
CNTRL2-GAIN*OVERO*COS ( SOS ) /DTR 
WRITE (105,61)  NREV , CNTRL1 , CNTRL2 , El , E2 
61  FORMAT (IX, 5G15 .8) 

GOTO  100 
ELSE 

OVERO-OVER 

ENDIF 

ELSE  IF(BEEP.EQ. 'Y' . AND. ABS(PSI -AINT(PSI/2 ,/PI)*2 . *PI-SOS) 

&. LT . DELPSI/1 . 9 .AND . (PSI-PSIO) ,GT.4.*PI)  THEN 
UNDER-BETA- BETLIM 

C  0.2  IN  THE  FOLLOWING  FORMULA  IS  THE  FLAPPING  ANGLE  ERROR  (IN  DEGREES)  ALLOWED 
IF ( (UNDER- .2*DTR) .GT.O. )  THEN 
U( 10) -ABS (BETA) 

DO  220  IU-1,9 
220  U(IU)-U0(IU+1) 

CALL  FDBACK(UNDER- . 2*DTR) 

CNTRL1- - GAI N* ( UNDER -0.2*DTR)*SIN(SOS) /DTR 
CNTRL2-GAIN* (UNDER- 0 . 2*DTR) *COS ( SOS ) /DTR 
WRITE (105,61)  NREV , CNTRL1 , CNTRL2 , El , E2 
SINGA1>'0' 

GOTO  100 

ENDIF 

ENDIF 

Ct  +  i  >  h  i  n  h  n  t  h  h  t  h  h  n  i  n  m  ntmi  h  mum 
C 

2790  PS I— PS 1+ DELPS I 

PSII-PSII+  DELPSI 
IF  (PSII . LE. 2 .*PI)  GO  TO  9801 
TAVG-TAVG/2 . /PI*B 
CTRW-TAVG/RHO/AREA/VT**2 
DO  9800  1-1,5 

CFAC-V**4+ ( TAVG /RHO /AREA ) ** 2+W** 3*2 . *V*SIN(Al+ALPHA) 

IF(CFAC.LT.O. )THEN 

W— SQRT(0 . 5*( -V**2+SQRT(V**4+(TAVG/RHO/AREA)**2) ) ) 

GO  TO  9800 
ENDIF 

WFAC— - V**2+SQRT ( CFAC ) 

IF(WFAC.LT.O. )THEN 

W-SQRT ( 0 . 5* ( - V**2+SQRT ( V**4+ ( TAVG/RHO/AREA) **2 ) ) ) 

GO  TO  9800 
ENDIF 

W-SQRT ( 0 . 5* ( - V**2+SQRT ( CFAC ) ) ) 

9800  WVT-W/VT 
TAVG— 0 . 

PSII-0. 

9801  IF( (PSI-PSIO) ,GT.NPRED*2 .*PI)  THEN 
BEEP-'N’ 

CNTRLl-0 . 
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CNTRL2-0 . 

E1-THETA1/DTR 

E2-THETA2/DTR 

WRITE <105, 61)  NREV, CNTRL1 , CNTRL2 , El , E2 
GOTO  100 
END  IF 
GOTO  2280 

100  WRITE(107 , 62)  NREV, GAIN 
62  FORMAT ( IX, 2G15 . 8) 

RETURN 

END 

C 

C - - - - - - - - - - - - - - - - - — - - - - - - - - . . 

C  This  subprogram  renews  the  GAIN  so  that  the  system  stability  may  be 
C  guaranteed  although  the  gain  so  obtained  makes  the  settling  time  longer. 
SUBROUTINE  FDBACK(OVER) 

REAL  U(10) ,U0(10) 

CHARACTER* 1  BEEP 

COMMON/PMT/PI , DTR , STEP , NMAX , PS 10 , DELPSI , W0 
COMMON/INPUT /THETA0 , THETAl , THETA2 , ALPHA , BETA0 , BETAD0 
COMMON/S IMU/B , D , CO , CT , EPS , KBETA , THETAT , MW , MIF , V , VT , DELX 
6. , RHO , KTEN , CRATEO , CRATE1 , CRATE2 , CRATEA , GAIN , SOS , BEEP , U , UO 
C 

C  An  adaptive  GAIN  is  chosen  so  a«5  to  guarantee  the  system  stability. 

C  The  design  principle  is  to  place  more  confidence  on  more  recent  controls 
C  with  exponential  confidence  coefficients. 

C  ANORM(U)  is  a  functional  subprogram  which  gives: 

C  ANORM(U)-[U(10)+exp(-c*l)*U(9)+. . .+exp(-c*9)*U(l)]/2-NORM  OF  U 
C  where  c  is  the  decaying  coefficient. 

C 

C  THE  FOLLOWING  IF-ENDIF  BLOCK  GIVES  A  SEMI -ADAPTIVE  GAIN 

IF (ANORM(U) .GE.ANORM(UO) . AND . ANORM(UO) .GT.l.E-4)  THEN 
GAIN-ANORM (U0 ) /ANORM (U ) *GAIN 
END  IF 

IF(GAIN.LE.O.S)  GAIN-0.8 
C 

THETAl-THETAl -GAIN*OVER*S IN ( SOS ) 

THETA2-THETA2+GAIN*OVER*COS ( SOS ) 

C  The  physical  limits  on  THETAl  and  THETA2  are  given  below: 

IF(THETAl . GE . 15 . *DTR)  THETAl-15 . *DTR 
IF (THETAl . LE . - 15 . *DTR)  THETAl— 15 . *DTR 
IF(THETA2.GE.15.*DTR)  THETA2-15 . *DTR 
IF(THETA2 . LE . - 15 . *DTR)  THETA2— 15 . *DTR 
DO  10  1-1,10 
10  U0(I)-U(I) 

RETURN 

END 

C 

C  i  -  - - -  - —  ...  -  - ■■  - 

C  This  subprogram  defines  a  special  kind  of  NORM  for  a  vector  used  in 
C  subprogram  FDBACK.  The  NORM  so  defined  has  the  following  properties: 

C  1.  NORM(U)  is  smaller  than  or  equal  to  inf inite-Norm  of  U 

C  2.  Each  coefficient  of  U  has  different  weighting  on  NORM(U), 

C  which  is  different  from  the  usual  P-Norm. 
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FUNCTION  ANORM(U) 

REAL  U(10) 

C-0.1 
ANORM-O . 0 

DO  ' D  1-1,10 

10  ANORM-ANORM+EXP ( - C*FLOAT ( 10 - I ) ) *U ( I ) 

RETURN 
END 
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Input  data  files  for  FORTRAN  programs 
(AH-1J  case) 
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FLP5.DAT  FOR  THE  AH-1J 

44. 

2.25 

2.25 

.01 


FLP6.DAT  FOR  THE  AH-1J 

0. 

-10. 

3122. 

1422. 

61.0 


FLP7.DAT  FOR  THE  AH-1J 


/  . 
15.27 
1.73 
0.11 
-4.48 
2.71 
-1.24 
2.6 


FLP8.DAT  FOR  THE  AH-1J 

.05 

5. 

.002378 

9500. 


CL. DAT  FOR  0012  AIRFOIL 


0. 

0. 

2. 

.211 

4. 

.422 

6. 

.633 

8. 

.844 

10. 

1.055 

11. 

1.161 

12. 

1.255 

13. 

1.334 

14. 

1.333 

15. 

1.19 

16. 

1.007 

21. 

.800 

39. 

1.18 

49. 

1.18 

129. 

1. 

147. 

1. 

161. 

.62 

172. 

.78 

180. 

.0 

CD. DAT  FOR  0012  AIRFOIL 


0. 

.008 

1. 

.0083 

2. 

.0085 

3. 

.0088 

4. 

.0093 

5. 

.01 

6. 

.011 

7. 

.0122 

8. 

.0138 

9. 

.0154 

10. 

.0174 

11. 

.0196 

12. 

.022 

13. 

.0264 

14. 

.038 

15. 

.102 

16. 

.155 

21. 

.332 

30. 

.562 

50. 

1.392 

60. 

1.66 

70. 

1.84 

80. 

1.96 

90. 

2.02 

100. 

2.02 

110. 

1.852 

120. 

1.652 

140. 

1.042 

160. 

.302 

165. 

.242 

170. 

.132 

175. 

.062 

180. 

.022 

